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SUMMARY 


\ 

The  method  of  lines  is  used  in  this  report  for  solving  one  linear, 
two  nonlinear  elliptic  boundary  value  problems  and  a  linear  eigenvalue 
problem.  An  analysis  of  the  stability  and  convergence  is  made  in  the  linear 
cases. 

1\ 


RESUME 


Dans  le  present  rapport,  on  exploite  la  methode  des  lignes  pour 
resoudre  un  probleme  lineaire,  deux  problemes  non  lineaires  de  valeurs 
elliptiques  aux  limites  et  egalement  un  probleme  de  valeur  propre  lineaire. 
On  effectue  une  analyse  de  la  stabilite  et  de  la  convergence  dans  les 
situations  lineaires. 
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APPLICATION  OF  THE  METHOD  OF  LINES  TO  THE  SOLUTION 
OF  ELLIPTIC  PARTIAL  DIFFERENTIAL  EQUATIONS 

CHAPTER  1.0  INTRODUCTION 

The  method  described  in  this  report  is  known  as  the  method  of  lines  ( from  here  on  we  refer 
to  the  method  as  MOL)  in  the  Soviet  Union,  where  it  has  been  used  for  some  forty  years.  The  basic 
feature  of  the  method  is  that  derivatives  with  respect  to  one  of  the  independent  variables  remain  con¬ 
tinuous,  while  derivatives  with  respect  to  the  other  independent  variables  are  replaced  by  finite- 
difference  approximations.  For  a  two-dimensional  problem  in  a  rectangle  the  region  could  be  consid¬ 
ered  as  divided  into  strips  by  dividing  lines  (hence  the  name)  parallel  to  one  of  the  axes.  At  each  line 
the  derivatives  normal  to  that  line  would  be  replaced  by  finite  differences  and  the  other  variable  left 
continuous.  Thus  the  system  of  partial  differential  equations  is  replaced  by  a  system  of  ordinary  dif¬ 
ferential  equations.  The  resulting  ordinary  differential  equations  may  then  be  solved,  at  least  in  some 
cases,  by  analytic  methods.  For  instance  Poisson’s  equation  with  linear  boundary  conditions  has 
received  much  attention,  Liskovets  (1965)  and  Leser  and  Harrison  (1966).  In  the  case  of  more  general 
equations,  particularly  those  of  nonlinear  type,  analytic  solutions  of  the  ordinary  differential  equations 
may  be  impossible  and  the  problem  must  be  treated  as  a  two-point  boundary-value  problem  to  be 
solved  numerically.  This  problem  may  then  be  solved  by  either  a  boundary -value  technique  such  as 
finite  differences  or  by  the  shooting  method  for  two-point  boundary-value  problems.  The  former 
technique  would  be  equivalent  to  solving  the  original  problem  by  the  grid  finite-difference  method. 

The  shooting  method  involves  estimating  unknown  conditions  at  the  initial  point  and  integrating  the 
ordinary  differential  equations  across  to  the  end  point.  The  required  boundary  conditions  at  the  end 
point  can  then  be  satisfied  by  iterating  on  the  missing  initial  conditions.  Because  of  the  elliptic  nature 
of  the  partial  differential  equations  this  initial-value  integration  is  strictly  improper.  Indeed  it  can  be 
shown  (Chapter  3)  that  the  ordinary  differential  equations  are  inherently  unstable.  One  of  the 
purposes  of  this  paper  is  to  convince  the  reader  that  in  many  physical  problems  of  interest  accurate 
solutions  can  readily  be  obtained  by  MOL  even  though  the  problem  is  incorrectly  posed.  It  is  shown 
that  if  the  region  of  interest  is  divided  into  sufficiently  few  strips  by  the  dividing  lines  then  accurate 
solutions  can  be  obtained  by  using  high-order  finite-difference  approximations.  As  more  and  more 
strips  are  taken  the  results  may  at  first  improve  but  they  will  eventually  become  meaningless  and  the 
iteration  technique  will  not  converge  to  a  solution. 

The  work  done  in  the  Soviet  Union  on  MOL  has  largely  been  limited  to  solving  linear 
equations  of  elliptic  (as  well  as  parabolic  and  hyperbolic)  type.  A  1965  review  paper  by  Liskovets 
gives  an  extensive  list  of  references  to  provide  the  mathematical  background  and  development  of  MOL. 
These  workers  have  developed  analytic  solutions  of  the  linear  ordinary  differential  equations  for 
certain  cases.  Also,  solutions  of  Poisson’s  equation  with  linear  boundary  conditions  were  obtained  in 
the  United  States  by  Leser  and  Harrison  (1966),  again  using  analytic  solutions  of  the  ordinary  differ¬ 
ential  equations. 

It  appears  that  MOL  (and  a  similar  technique  called  the  method  of  integral  relations)  was 
first  used  in  nonlinear  problems  for  the  supersonic  blunt -body  problem  which  is  of  interest  to  aero- 
dynamicists,  Belotserkovskii  (1965).  Klunker,  South  and  Davis  (1971)  have  discussed  more  recent 
applications  of  the  method  to  the  solution  of  equations  of  elliptic  type  such  as  the  supersonic  blunt- 
body  problem  and  conical  flow  problems  which  are  of  great  importance  in  aerodynamics.  In  general 
the  method  has  received  more  attention  for  solving  the  correctly  posed  parabolic  type  of  equation. 
Aktas  (1978)  gives  a  recent  review  of  some  applications  of  MOL  to  parabolic  and  hyperbolic  as  well  as 
elliptic  problems. 

After  describing  the  method  of  lines  in  Chapter  2  we  then  carry  out  a  stability  analysis  in 
Chapter  3.  Some  examples  are  next  discussed  in  Chapter  4  and,  finally,  in  Chapter  5,  we  look  at 
solutions  to  eigenvalue  problems. 


-2- 


REFERENCES  FOR  CHAPTER  1 


Aktas,  Z. 


Belotserkovskii,  O.M. 
Chushkin,  P.I. 


Klunker,  E.B. 
South,  J.C.,  Jr. 
Davis,  Ruby  M. 

Leser,  T. 
Harrison,  J.T. 


Liskovets,  O.A. 


On  the  Application  of  the  Method  of  Lines. 

Second  International  Conference  on  App.  Num.  Modelling, 

Sept.  1978. 

The  Numerical  Solution  of  Problems  in  Gas  Dynamics. 

IN  Basic  Developments  in  Fluid  Dynamics  (M.  Holt,  Ed.),  Vol.  1, 
Academic  Press,  New  York,  1965. 

Calculation  of  Nonlinear  Conical  Flows  by  the  Method  of  Lines. 
NASA  TRR-374, 1971. 


The  Method  of  Lines  for  Numerical  Solution  of  Partial 
Differential  Equations. 

Ballistic  Research  Laboratories,  BRLR  1311,  March,  1966. 

The  Method  of  Lines  (Review),  Differential  Equations  1(1 965). 
Translation  of  Differential’nye  Uravneniya  1  (1965),  pp.  1308- 
1323. 


CHAPTER  2.0  THE  METHOD  OF  LINES 


FIG.  2.1  MOL  APPLIED  TO  A  RECTANGLE 


This  simple  linear  case  is  sufficiently  detailed  to  describe  the  principles  of  the  method  of  lines.  More 
complicated  equations  and  boundary  conditions  can  be  solved  in  the  same  manner.  Minor  differences 
in  the  treatment  of  other  equations  and  boundary  conditions  are  seen  in  the  examples  of  Chapter  4. 
For  example,  mixed  derivatives  and  higher  order  partial  differential  equations  are  handled  in  the 
Example  of  Section  4.2. 


Other  geometries  may  be  treated  after  transformation  as  shown  later. 

We  first  take  equally  spaced  lines  parallel  to  the  x  axis  and  number  them  (see  Fig.  2.1) 

b 

0,1,2, .  .  .  N  where  N  is  the  number  of  divisions  in  the  region  0<y<b,  thus  h  =  — . 

N 

Next  we  write  the  Equation  (2.1)  as  a  set  of  ordinary  differential  equations.  Letting  p  ~  4>x 
we  have,  at  each  line  i  (i  =  0,1  .  .  .  N- 1) 


—  =  f(x,yf) 

dx 


'J'.+i  -  Wi  +  *i-\ 


(2.4a) 


(2.4b) 


since  we  have  approximated  i//vv  by  the  three  point  difference  formula 


^  _  <Hx,y+h)  -  2i//(x,y )  +  ^(x,y-h)  t  0(h2^„,) 

h2 


At  line  N  we  have 


i//N  =  sinux 


(2.6a) 


while  at  the  image  line  i  =  -1  we  have 


(2.6b) 


The  initial  and  end  conditions  for  the  system  are 


^  =  0 


at  x  =  0  and 


Pi  =  0 


at  x  =  a  to  ensure  symmetry  of  \p  about  that  point. 

The  system  given  by  (2.4)  is  now  a  system  of  ordinary  differential  equations  with  two  point 
boundary  conditions.  Thus  it  can  be  solved  by  standard  shooting  techniques  used  for  solving  two 
point  boundary  value  problems.  Some  of  these  methods  are  given  in  Keller  (1968)  in  which  the 


f-v 
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contracting  map  approach  and  the  variational  approach  using  Newton’s  method  of  iteration  are 
described.  The  latter  method  is  much  faster  than  the  former  and  so  is  to  be  recommended.  We 
describe  next  methods  which  are  similar  to  the  variational  approach  but  which  do  not  set  up  the 
variational  equations  as  such.  They  are  thus  easier  to  program  and  yet  still  have  the  advantage  of 
second  order  convergence.  We  first  discuss  Newton’s  method  and  then  Powell’s  (1965)  method  as 
applied  to  the  shooting  technique. 

2.2  Newton’s  Method  of  Iteration  Applied  to  the  Two  Point  Boundary  Value  Problem 

We  first  notice  that  the  problem  is  solved  completely  provided  that  we  know  the  initial 

slopes 


Pi 


=  Fs 


(i  =  0,1  .  .  .  N-l)  say  at  x  =  0.  This  is  clear  since  a  knowledge  of  F(  enables  us  to  integrate  the 
Equations  (2.4)  from  x  =  0  to  x  =  a  by  standard  methods  such  as  Runge  Kutta.  The  present  authors 
normally  use  Hamming’s  predictor  modifier  corrector  (PMC)  method,  see  Hamming  1959,  with  the 
Runge  Kutta  starting  procedure.  The  PMC  method  has  a  discretization  error  0(5xs )  and  requires  half 
as  many  gradient  evaluations  as  the  full  Runge  Kutta  method.  On  each  evaluation  of  the  gradients  in 
(2.4)  the  boundary  conditions  (2.6)  are  used  and  thus  automatically  satisfied.  Hence  integration  to 
x  =  a  (provided  round  off  errors,  discretization  errors  and  inherent  instability  are  negligible  —  see  later) 
is  achieved  and  we  recover  the  boundary  condition  (2.8).  Now,  all  boundary  conditions  are  satisfied 
as  well  as  the  differential  equations  and  we  have  a  complete  solution. 

However,  since  we  do  not  know  F;  a  priori,  we  must  develop  some  scheme  for  improving  on 
a  given  estimate  of  Fj.  We  notice  that  the  only  boundary  conditions  not  satisfied  after  an  integration 
of  (2.4)  with  a  given  estimate  of  Fj  are  the  conditions  (2.8).  The  procedure  to  follow  therefore  is  to 
di/'j 

drive  the  values - at  x  =  a  closer  to  zero  by  suitably  adjusting  Fj. 

rlv  * 


Let 


d^j 

dx 


(2.9) 


at  x  =  a.  Then  hef  can  be  made  smaller  by  making  changes  to  Fj  as  indicated  by  the  Newton  iteration 


N-l  de; 

2  I^r5Fj  =  =  0,1...  N-l)  (2.10) 

j=0  d*i 


The  above  form  requires  approximations  to  the  partial  derivatives 


3ej 

3F, 


To  obtain  these  Keller  (1968),  amongst  others,  uses  the  variational  equations.  However  the  present 
authors  have  found  that  it  is  not  necessary  to  use  the  variational  equations,  instead  one  can  write 
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ei(FoF!  .  .  .  Fj  +  AFj  . . .  FN_j)  -  ej(F0Fj  . . .  Fj  .  .  .  Fn_j) 
3F^  ~  AFj 


(2.11) 


where  AF,  is  a  small  value  normally  taken  as  10'6  AFj  or  10-6  if  | AFj | <  1  (64  bit  words  with  approxi¬ 
mately  16  decimal  digit  accuracy  have  been  used  in  all  our  computations).  In  (2.11)  the  6j  have  been 
written  as  functions  of  F0  .  .  .  Fn_j  since  a  knowledge  of  F0  . .  .  Fn_j  enables  us  to  integrate  (2.4) 
from  x  =  0  to  x  =  a  and  so  we  find  the  as  implicit  functions  of  F0  .  .  .  Fn_j. 

Notice  that  the  form  (2.11)  is  exact  if  6;  is  a  linear  function  of  F0  .  . .  Fn_,  . 

The  procedure  to  find  the  derivatives 


3ej 


is  therefore  to  make  an  integration  of  (2.4)  with  the  latest  estimate  F0  .  .  .  Fn_j  .  Then  make  successive 
integrations  with  Fj  changed  to  F=  +  AFj  for  j  =  0,1  . . .  N-l  and  after  each  integration  substitute  into 
(2.11)  to  find  the  required  partial  derivatives.  Having  completed  these  integrations  we  now  substitute 
(2.11)  into  (2.10)  and  find  the  changes  6Fj  by  standard  linear  methods  such  as  Gaussian  elimination. 

The  values  thus  obtained  for  6Fj  can  now  be  used  to  improve  on  the  last  estimate  for  Fj. 

This  use  of  5Fj  is  described  below. 

2.2.1  The  Linear  Case 

If  the  differential  equations  and  boundary  conditions  are  linear  as  in  our  example  then  new 
values  of  Fj  given  by 


Fj  =  Fj  (old)  +  5Fj 


with  SFj  found  from  (2.10)  are  the  correct  values.  This  follows  from  the  fact  that  Newton’s  scheme 
(2.10)  is  quadratic  in  convergence  so  that  if  fj  is  a  linear  function  of  Fj  then  the  correction  6Fj  is 
exact.  In  addition  the  formula  (2.11)  is  exact  for  linear  systems. 

2.2.2  The  Nonlinear  Case  and  the  Modified  Newton  Method 

It  is  well  known  that  Newton’s  method  will  not,  in  general,  converge  if  the  functions  e,  are 
nonlinear  functions  of  F:  unless  the  initial  estimate  is  sufficiently  near  the  solution.  However  it  can  be 
shown  (see  for  example  kowalik  and  Osborne,  1968)  that  the  direction  given  by  (2.10)  i.e. 


6F  =  (6F0,  5F,, . . .  5Fn.,) 

j 

is  downhill  in  the  sense  that  Ze?  will  decrease  by  changing  Fj  such  that 

Fj  =  Fj  (old)  +  XSFj  (2.12) 

with  X>0.  This  feature  leads  us  to  the  modified  Newton  procedure  whereby  X  in  (2.12)  is  chosen  such 
that  Zef  using  new  values  of  Fj  is  less  than  the  previous  Zef  using  old  values.  Thus  we  select  X  =  1, 
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•y 

compute  Fj  from  (2.12)  and  hence  compute  new  values  of  ZCj  by  integration  as  described  previously 
and  accept  the  new  values  of  Fj  if  they  give  an  improvement  to  £e?.  If  the  new  values  do  not  yield  an 
improvement  we  take  X  =  1/2  and  repeat  the  procedure.  On  successive  failures,  if  they  occur,  we  use 
X  =  1/4, 1/8  .  .  .  and  so  on  until  success  is  achieved.  Once  we  have  a  success  we  then  use  (2.10)  to  find 
a  new  direction  in  which  to  advance.  This  method  has  been  used  successfully  by  Jones  (1973),  amongst 
others. 


Some  improvement  in  efficiency  can  be  made  to  the  above  modified  Newton  method  and 
this  is  now  described.  Newton’s  method  requires  a  knowledge  of  the  partial  derivatives 


3ej 

and  this  determination,  from  (2.11),  may  be  costly  since  N  integrations  of  (2.4)  are  needed.  In 
Newton’s  method  as  described  above  these  partial  derivatives  have  to  be  evaluated  at  each  step  given  by 
(2.10).  However  Newton’s  method  can  in  some  cases  be  improved  on  by  using  the  following  simplified 
version  (see  Collatz  (1966)  for  more  details). 

Essentially  the  method  consists  of  using  (2.10)  repeatedly  with  values 


3F: 


unchanged  from  an  earlier  step.  The  simplified  version  is  illustrated  in  Figure  2.2  for  the  one  unknown 
F0.  Clearly 


FIG.  2.2  SIMPLIFIED  NEWTON  METHOD 


convergence  measured  in  terms  of  number  of  steps  is  slower.  However  it  has  the  advantage  that  the 
derivatives 


do  not  have  to  be  evaluated  repeatedly  nor  does  the  system  (2.10)  have  to  be  solved  completely  since 
the  inverse  of 


J  = 


9e 

9F 


can  be  computed  once  and  stored.  South  and  Klunker  (1969)  have  used  this  method  successfully  for 
conical  flow  calculations.  They  find  that  the  method  works  well  if  one  is  sufficiently  near  the  solution. 
If  not  one  may  use  several  modified  Newton  steps  prior  to  using  the  simplified  version. 

2.3  Powell’s  Method  of  Iteration 

Powell’s  method  for  minimizing  a  sum  of  squares  of  nonlinear  functions  is  given  in  Powell 
(1965).  The  method  minimizes 


N-l 

2 

i=0 


with  respect  to  F0)F,  .  .  .  Fm_,  (M<N),  where  the  N  functions  e,  are  nonlinear  functions  of  the  M 
unknowns  F: .  The  method  is  essentially  that  of  least  squares  minimization  in  which  Ze*  is  minimized 
by  making  changes  to  Fj  indicated  by  the  direction  6F  given  by 


M  l  (N-l 

9ek  9ek 

N-l  9ek 

2 

{  2 

•  SFj  =  -  2  ek  —  ( 
k=0  aF> 

i  =  0  .  . 

.  .M-l) 

(2.13) 

3=0 

[k=0 

9Fj  9Fj 

The  step  to  find  new  values  of  F  =  (Fg^j  .  .  .  FM-1 )  is  given  by 


F  =  F  (old)  +  X5F 


in  which  X  is  chosen  (by  search)  such  that  Ze?  is  minimized  along  the  direction  6F.  During  the  search 
along  6F  to  locate  the  minimum,  functions  e;  have  to  be  evaluated  at  different  values  of  X;  thus  one 
can  calculate  an  approximation,  by  differences,  to  the  rate  of  change  of  Cj  along  the  direction  SF  at 
the  new  minimum  point.  Powell  shows  how  these  partial  derivatives  can  be  used  in  conjunction  with 
previous  values  of 


9ek 

9?; 


to  determine  the  next  step  given  by  (2.13). 
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In  principle  the  method  guarantees  convergence  since  we  do  not  take  a  step  unless  Sej2 
decreases.  In  practice  the  method  has  been  found  to  fail  but  only  rarely  and  usually  when  a  poor 
estimate  has  been  made.  Powell’s  method  when  it  first  appeared  in  1965  was  probably  the  most 
efficient  method  of  its  day.  Other  methods  which  may  now  be  more  efficient  than  Powell’s  are  the 
spiral  method  of  Jones  (1970)  and  Peckham’s  (1970)  method  but  these  have  not  been  tried  on  two 
point  boundary  value  problems. 

The  advantage  of  Powell’s  method  over  the  modified  and/or  simplified  Newton  method  is 
firstly  that  Powell’s  method  is  usually  more  efficient  and  secondly  that  fewer  “unknowns”  can  be  used 
to  determine  the  solution.  In  MOL  it  may  be  that  the  unknown  function,  in  our  example 


at  x  =  0,  can  be  represented  in  a  series  expansion,  say 


—  =  S  Fjcosjy 


(M«N)  and  in  this  way  the  work  involved  in  finding  the  partial  derivatives 


(2.14) 


is  greatly  reduced  since  now  we  have  to  integrate  (2.4)  only  M  times  to  obtain  the  required  partial 
derivatives.  Such  a  form  (2.14)  was  used  in  the  conical  flow  calculation  by  Jones  (1968). 

Powell’s  method  also  has  quadratic  convergence  provided  one  is  sufficiently  near  the 
solution  and  provided  e;  =  0  at  the  minimum.  Note  that  the  correct  solution  is  again  obtained  in  one 
step  if  the  system  is  linear. 

During  the  preceding  discussion  the  phrase  “provided  one  is  sufficiently  close  to  the 
solution”  keeps  recurring  and  indeed  it  is  very  important  in  nonlinear  cases  to  have  reasonable 
estimates  of  the  unknowns  particularly  to  cut  down  computing  costs.  In  view  of  the  importance  we 
discuss  initial  estimates  in  a  separate  section  to  follow. 

2.4  Initial  Estimates 

2.4.1  The  Linear  Case 

In  the  case  where  the  differential  equations  and  boundary  conditions  are  linear  the  MOL 
equations  can  be  solved  in  one  step  using  Newton’s  or  Powell’s  method  whatever  initial  estimate  is 
made. 

However  the  above  statement  must  be  viewed  with  some  caution.  Although  it  is  true  in 
principle,  in  practice  we  clearly  have  to  have  an  initial  estimate  which  at  least  allows  integration  to  the 
end  boundary  x  =  a  without  the  solution  blowing  up.  In  the  authors’  experience  it  is  possible  to  take 
quite  crude  estimates  such  as,  at  x  =  0, 
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_  ffy 
dx  b 


(2.15) 


in  the  example  of  Section  3.5.  Compare  this  to  the  exact  solution 


d^  coshiry 

-  =  7T - 

dx  coshffb 


The  motivation  behind  the  choice  (2.15)  was  merely  that  it  satisfied  the  boundary  condition  at  y  =  b 
evaluated  at  x  =  0. 

2.4.2  The  Nonlinear  Case 

As  has  been  pointed  out  it  is  desirable  (and  sometimes  essential)  to  have  a  reasonably  good 
estimate  in  the  nonlinear  case.  The  authors  have  found  that  this  limitation  is  not  severe  for  two 
reasons. 


The  first  reason  is  that  a  nonlinear  problem  can  often  be  made  linear  by  a  suitable  choice  of 
a  parameter  in  the  problem.  This  linear  problem  can  then  be  solved  with  a  fairly  crude  estimate  and 
then  the  parameter  can  be  varied  in  discrete  steps.  To  obtain  the  solution  at  each  value  of  the  para¬ 
meter  a  good  initial  estimate  is  available  by  extrapolation  from  previous  results. 

The  second  reason  is  that  a  parameter  in  the  problem  can  often  be  chosen  such  that  a 
solution  is  already  known  at  that  value  and  estimates  for  each  successive  values  of  the  parameter  are 
then  obtained  by  extrapolation  as  above.  An  example  of  this  is  in  first  setting  the  angle  of  incidence 
to  be  zero  in  the  conical  flow  calculations  of  Section  4.3;  this  has  the  effect,  for  the  circular  cone,  of 
making  the  flow  axisymmetric  and  solutions  in  this  case  are  well  known. 

Note  that  such  an  extrapolation  (the  authors  use  quadratic  extrapolation  as  soon  as  three 
previous  solutions  are  available)  is  not  restricted  by  computer  storage  limitations  as  it  may  be  in  grid 
techniques.  This  follows  since  we  have  only  0(N)  unknowns  to  store  rather  than  0(N2). 

2.5  Termination  Procedures 

Good  termination  procedures  for  nonlinear  problems  are  often  difficult  to  find  and  some 
attention  should  be  paid  to  them.  Clearly  in  the  linear  case  one  step  is  all  that  is  needed  and  iteration 
can  then  cease.  Some  of  the  subprograms  written  for  Powell’s  method  terminate  when  the  next 
change  to  Fj,  5Fj,  is  less  than  a  certain  amount  say  [6Fj|<p  for  some  small  value  of  p.  However  in 
some  cases  a  better  criterion  might  clearly  be 


5Fi 


<  P 


There  may  be  a  problem  in  choosing  p  since  too  small  a  value  may  result  in  excessive  computer  time. 


The  authors  have  found  it  better  to  use  a  criterion  which  depends  on  the  size  of  the  residuals 
of  6j  —  either 


1_ 

N 


N-l 

1  <  5 

0 
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or 


max|ej|<5 

i 


for  example  can  be  used.  In  the  case  of  the  conical  flow  calculations  of  Section  4.3  the  choice  of 
residual  could  be  found  from  physical  reasoning  and  was  related  to  the  mass  flow  into  or  from  the 
body  relative  to  the  total  free  stream  mass  flow  being  less  than  a  certain  amount. 

A  combination  of  the  above  criteria  together  with  a  preset  maximum  number  of  integrations 
is  normally  employed. 

2.6  Other  Systems 

In  order  to  describe  MOL  a  simple  linear  example  was  discussed.  However,  the  procedure  of 
solution  for  other  systems  of  equations  is  identical  to  that  already  described,  particularly  as  we  have 
already  covered  Newton’s  and  Powell’s  methods  for  solving  nonlinear  problems.  The  user  may  need  to 
use  other  finite  difference  formulas,  for  example  for  a  first  derivative  we  could  write 


di// 

3y 


»(*,y+h)  -  »(>,y-h)  +  (Hh^1,I) 


(2.16) 


Formulas  (2.5)  and  (2.16),  accurate  to  0(h2),  which  use  values  on  3  adjacent  lines,  may  be 
used  in  MOL  applications.  However  a  decided  improvement  can  be  obtained  by  using  formulas 
obtained  from  values  on  five  adjacent  lines.  These  formulas  are 


d\J/  4  \p(x, y+h)  -  Mx,y-h) 

6y  3  2h 


1  Mx,y+2h)  -  ^(x,y-2h) 
3  4h 


+  0(h4tf/v) 


(2.17a) 


and 


=  £  Mx,y+h)  +  ^(x,y-h)  -  2\p(x,y) 
dy2  3  h2 


1  ift(x,y+2h)  +  Mx,y-2h)  -  2yfr(x,y) 
3  4h2 


+  0(h4^VI) 


(2.17b) 


The  advantages  of  using  the  five  line  schemes  (2.17)  in  place  of  the  three  line  schemes  is 
shown  below.  The  next  section  then  gives  difference  formulas  of  the  same  accuracy  0(h4)  as  (2.17) 
which  can  be  used  on  the  boundary  or  at  the  line  adjacent  to  the  boundary. 


-  -iisrv 
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Suppose  E  represents  the  exact  derivative  of  a  function  \p  with  respect  to  y  and  let  Aj  and  A2 
represent  the  approximations  given  by  Equations  (2.16)  and  (2.17a)  respectively.  Then  it  can  be 
shown,  by  looking  at  the  truncation  errors,  that 


A1 

8y,2 

*1  = 

1  — 

—  - 

E 

6 

and  that 


t<  A2  6y24 

62  =  I1"  ¥  ~~2A  ^7 


where  5y ,  and  8 y2  are  the  finite-difference  increments.  Now  we  want  to  find  the  ratio  6y2 :6y  2  to  give 
the  same  accuracy  in  both  formulas,  i.e.,  =  e2  =  e.  For  this  condition  we  have 


5y2  ^iii2  - 

—  -  0.90  — -  e~ 


(2.18) 


The  value  of  e’/48y2/8y1  can  be  calculated  from  this  formula  for  well-known  functions.  Its 
value  is  approximately  0.90  for  sin  ny,  cos  ny  and  exp(ny),  while  for  log  y  its  value  is  0.57.  Hence 
provided  the  approximate  formula  (2.18)  holds,  it  can  be  seen  that  if  e  is,  say  104  (0.01%  accuracy), 
then  8y2/8y!  lies  between  6  for  the  log  function  and  9  for  the  sin,  cos  and  exp  functions. 

The  above  analysis  shows  that  6-9  times  as  many  dividing  lines  must  be  used  with  (2.16)  to 
get  equivalent  accuracy  to  (2.17a)  for  the  first  derivatives.  Equation  (2.16)  requires  about  half  as 
many  computer  multiplications,  divisions,  etc.,  compared  to  (2.17a)  but  this  affects  only  one  state¬ 
ment  of  the  computer  program  and  so  is  insignificant  in  terms  of  computer  time.  A  similar  saving  in 
lines  may  be  made  by  using  the  5-point  scheme  (2.17b)  instead  of  (2.5)  for  second  derivatives.  The 
superiority  of  the  five  line  schemes  over  the  three  line  schemes  will  be  illustrated  in  Section  3.5. 

2.7  Difference  Formulas  Applied  Near  a  Boundary 

To  apply  (2.17)  on  a  line  adjacent  to  a  boundary  line  is  not  possible  unless  the  boundary 
line  is  a  line  of  symmetry  when  image  lines  are  used.  If  the  boundary  line  is  not  a  line  of  symmetry 
then  the  following  formulas  are  recommended. 

(i)  Dirichlet  Boundary  Condition 


1  3  1 


-  ~ ('J'l-'M  +  0(h5  \j/v ) 


(2.19) 
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for  a  first  derivative,  where  lines  0,  1,  2,  3,  4  are  adjacent  and  line  0  forms  the  boundary.  The  symbols 
\p0,  \p  |,  etc.,  refer  to  the  values  of  \p  on  line  0,  line  1,  etc.,  and  \pv  refers  to  the  fifth  derivative  of  \J/ 
with  respect  to  y  in  the  range  considered.  The  upper  sign  is  used  if  the  lines  0, 1,  2,  3,  4  are  at  increas¬ 
ing  y  values,  otherwise  the  lower  sign  is  used.  For  a  second  derivative 


4 


+  0(h6^vi)  (2.20) 


h2  3  Vi  l  1.75 

7> — 7  =  -  t(^2-^i)  +  -r~(^3-^i) 
2  3  2  6  3 


0.125  ,  ,  1.25  , 

+  —7— (^5-iM  +  —T-ib-ti) 


is  recommended. 


(ii)  Neumann  and  Mixed  Boundary  Conditions 


In  this  case  9i///9y  is  a  constant  on  the  boundary  or  else  9i///9y  is  given  as  a  function  of  \j/. 
In  the  latter  case  calculate  d\p0/dy  from  the  boundary  condition  and  then,  for  both  cases,  use 


3^i  8.5  1  1  1  Wo  ,  v 

±h_r“  =  -  Ti'l'i-'l'i)  +  7^^r^3>  +  +  °(h  'Z'  )  (2.21) 


9y 


18 


3  9y 


for  a  first  derivative  and 


h2  3  3.5  1  0.0625 

2  dy2  4  9  6 


8.03125  0.625  3^o  ,  V1 

+  - Wo-ti)  ±  1  —  h  —  +  0(h6^vl) 

a  oy 


9 


(2.22) 


for  a  second  derivative.  The  appropriate  ±  sign  is  chosen  in  the  same  way  as  in  (i)  above. 
Also  the  appropriate  form  for  92i//0/9y2  is  given  by 


h2  32  *0  3^o  Wo 

—  — -  =  4 (^-^o  +  h— )  -  1.5(i^2-^o  +  2h— ) 

2  dy2  9y  3y 

4  1  3^o 

+  *  3h— )  -  *0  ?  4h— ) 


+  0(h6*VI) 


(2.23) 
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In  addition  to  first  and  second  derivatives,  formulas  for  higher  derivatives  may  be  found  if 
required.  For  example,  in  the  case  of  the  Neumann  boundary  condition,  93 \i/  j  /9y3  can  be  determined 
by  solving  the  equations 


(nh)2  _  (nh)3  .  (nh)4 


*n+i  -  =  nh^J  + 


h9^o  ,  , 

—  =  h*  -  h2* 

3y  1 


ii 

i 


h3  m 
+  —  \l/  11 

o  ’'l 


for  ^J11  in  terms  of  \J/0,  \jj j,  \p2,  ^3,  ^4  and  30o/3y.  Note  that  the  required  formula  can  be  found 
conveniently  by  computer  matrix  inversion. 


2.8  Other  Geometries 

For  simplicity  we  used  a  rectangle  to  describe  MOL.  The  application  to  the  geometry  of 
Figure  2.3  with  polar  co-ordinates  is  obvious;  in  this  case  lines  </>  =  constant  are  used. 


FIG.  2.3  MOL  WITH  POLAR  CO-ORDINATES 


Application  of  the  method  to  more  complicated  geometries  is  generally  limited  only  by  the 
degree  of  ingenuity  of  the  user,  since  an  appropriate  choice  of  co-ordinates  and  various  transformations 
can  be  used  to  map  most  regions  into  one  of  the  simple  regions  just  mentioned.  For  example  the 
region  ABCD  of  Figure  2.4  can  be  mapped  into  the  rectangle  0<£<1, 0<r?<l  by  the  transformation 


x-gj(y)  y-fi(x) 

g2(y)-gi(y)’  17  f2(x)-fi(x) 


4 


y 


=  o2( y ) 


c 


FIG.  2.4  THE  REGION  ABCD 


This  transformation  is  now  applied  to  the  partial  differential  equations  using,  for  example, 


and 


d_ 

9x 


d_ 

dr) 


d2 


a2 


d2 


—  =  Zxx  7Z  +  Vxx  T  +  K  +  <  ~  +  2^X  Vx 


dx2 


3* 


dr) 


0f2 


dr)1 


d2 

dtfr) 


and  then  the  resulting  equations  are  solved  by  MOL  in  the  (£,77)  plane.  Note  that  the  transfor¬ 
mation  preserves  linearity.  Such  a  transformation  is  used  to  solve  the  conical  flow  problem  of 
Section  4.3. 


2.9  A  Cautionary  Note  on  the  Instability 


The  preceding  description  of  MOL  is  somewhat  idealized  in  the  sense  that  we  have  assumed 
that  the  computer  round  off  errors,  the  discretisation  error  in  the  PMC  integration  scheme  and  also 
the  instability  inherent  in  the  ordinary  differential  equations  are  negligible.  They  of  course  are  not 
negligible;  in  fact  it  is  known  that  MOL  solutions  are  unstable  since  as  more  and  more  lines  are  used 
the  solutions  diverge  and  thus  become  meaningless.  This  fact  may  deter  people  from  using  MOL  and 
one  of  our  main  purposes  in  this  report  is  to  show  people  that  problems  of  physical  interest  have  been 
solved  successfully  by  MOL.  We  also  present  an  analysis  of  the  instability  in  the  linear  case  in  the  next 
chapter.  This  analysis  will  show  the  form  of  the  instability  and  thus  the  user  will  be  made  aware  of 
what  problems  to  expect.  Then  the  same  linear  problem  that  is  analysed  is  solved  numerically  by 
MOL.  It  will  be  seen  that  4  or  5  significant  figure  accuracy  is  obtained. 

Perhaps  we  should  point  out  for  the  moment  that  MOL  solutions  may  be  compared  to 
asymptotic  series  solutions  since  the  latter  are  divergent  series  so  that  as  more  and  more  terms  are 
taken  the  solutions  become  meaningless.  However  many  useful  results  are  obtained  with  asymptotic 
series  using  only  a  few  terms. 

Also  parallel  shooting  methods  (Keller  1968)  may  be  used  to  overcome  the  instability.  In 
these  methods  we  subdivide  the  total  range  of  integration  0<x<a  into  subintervals  and  ‘shoot’,  for 
example,  as  shown  diagrammatically  in  Figure  2.5  below.  Of  course  the  number  of  unknowns  is 
increased  in  this  case  since  we  have  to  estimate  values  at  x  =  0,  x  =  a/2  and  x  =  a.  In  the  linear  case  the 
total  computational  time  will  be  practically  the  same  as  using  the  one  complete  range  since  the  number 

of  unknowns  is  increased  by  a  factor  of  four  but  each  integration  to  get  the  derivatives- —  is  over  the 

dFj 

range  a/4.  Of  course  the  time  involved  in  solving  the  linear  system  (2.10)  or  (2.13)  would  eventually 
become  large  if  many  subintervals  were  used. 

The  parallel  shooting  method  is  used  on  our  linear  example  in  Section  3.5.2. 


H - 1 - 1 - 1 - h 

0  a/2  a 


x 


FIG.  2.5  PARALLEL  SHOOTING  RANGES 
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CHAPTER  3.0  STABILITY  AND  CONVERGENCE  OF  MOL 
3.1  Introduction 

In  the  first  part  of  this  chapter  Hadamard’s  example  of  the  Cauchy  problem  for  the  Laplace 
equation  in  a  rectangle  is  discussed.  This  example  shows  that  we  may  expect  MOL  to  be  inherently 
unstable. 


A  closer  analysis  of  MOL,  using  the  scheme  (2.4),  for  the  Laplace  equation  in  a  rectangle  is 
then  made.  This  analysis  shows  that  the  general  solution  by  MOL  is  comprised  of  two  parts.  The  first 
of  these  is  an  unwanted  solution  which  is  negligible  for  sufficiently  small  x  (the  continuous  variable)  or 
when  sufficiently  few  lines  are  used  but  which  otherwise  grows  large.  The  second  part  of  the  general 
MOL  solution  is  the  required  solution  which  tends  to  the  exact  solution  as  the  y  increment  tends  to 
zero  (i.e.  as  the  number  of  lines  increases). 

In  the  final  part  the  stability  and  convergence  of  MOL  is  illustrated  with  a  linear  example. 

It  will  be  seen  that  good  accuracy  (4  significant  figures)  is  obtained  with  a  suitable  choice  of  both 
number  of  lines  and  integration  step  size  5x. 

3.2  Hadamard’s  Example 

Hadamard  (1952)  investigates  the  solution  of  Laplace’s  equation  and  shows  that,  subject  to 
Cauchy  data  of  a  certain  type,  the  solution  is  not  well  behaved  since  it  will  oscillate  between  very  large 
positive  and  negative  values  when  the  correct  solution,  in  the  limit  of  vanishing  Cauchy  data,  should  be 
zero.  Hadamard  poses  the  example 


Mx  +  ^yy  =  0 


(3.1) 


with  Cauchy  data  given  at  the  line  x  =  0, 

MO,  y)  =  0, 

lM0,y)  =  An  sin  ny,  (3.2) 

where  n  is  large  and  An  is  a  function  of  n  which  grows  small  as  n  grows  large  (e.g.,  n  p,  p  >  0).  The 
solution  to  this  problem  is 

Mx,y)  =  (A„/n)  sin  ny  sinh  nx  (3.3) 


The  sinh  nx  factor  is  large  because  of  the  growth  of  enx.  The  sin  ny  factor  causes  oscillation  of  the 
function  with  varying  y.  Hence  however  close  to  zero  we  choose  to  make  the  Cauchy  data  (i.e.,  n  large) 
the  solution  i|/(x, y)  will  not  be  zero  but  will  oscillate  between  large  positive  and  negative  values.  Since 
zero  is  the  solution  of  (3.1)  with  vanishing  Cauchy  data  ( An  =  0)  we  conclude  that  for  Laplace’s 
equation  the  dependence  of  the  solution  on  the  initial  data  is  not  in  general  continuous. 

Garabedian  (1964)  concludes  also  that  the  above  problem  is  not  correctly  set  or  well  posed. 
He  defines  a  boundary-value  problem  for  a  partial  differential  equation,  or  for  a  system  of  partial 
differential  equations,  to  be  correctly  set  in  the  sense  of  Hadamard  if  and  only  if  its  solution  exists,  is 
unique,  and  depends  continuously  on  the  data  assigned. 
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Consider  now  Hadamard’s  example  in  the  context  of  MOL.  We  may  proceed  with  MOL  by 
estimating  <//x(0,y)  and  using  this  estimate  integrate  (3.1)  numerically  away  from  x  =  0.  After  iteration 
we  arrive  at  a  numerical  solution  of  i/'x(0,y).  Suppose  the  exact  solution  to  the  problem  is  zero,  i.e., 
i//x(0,y)  =  0,  but  that  due  to  discretization  and  round-off  errors  the  value  of  i//x(0,y)  is  of  order  10” 10. 
Then  the  situation  is  similar  to  n  being  large  in  (3.2)  in  the  sense  that  An  is  not  quite  zero.  As  we 
integrate  numerically  in  MOL  away  from  x  =  0  the  solution  will  be  in  a  form  similar  to  (3.3)  and  will 
thus  become  large  and  oscillatory  for  sufficiently  large  x.  Thus  we  cannot  obtain  a  good  approximation 
to  the  exact  solution  unless  x  is  small. 

It  may  also  be  noted  that  in  the  general  case  when  solving  by  MOL  even  an  exact  i//x(0,y) 
cannot  give  a  good  solution  for  all  x.  The  reason  for  this  is  that  the  x  discretization  error  introduces 
an  unwanted  solution  equivalent  to  An  f  0  in  the  Hadamard  example.  Thus  the  instability  is  always 
present  but  its  contribution  may  be  insignificant  for  x  sufficiently  small. 

The  above  observations  of  instability  are  analyzed  more  closely  in  the  next  part  of  this 
chapter.  It  will  be  confirmed  that  a  reasonably  accurate  solution  may  be  obtained  if  x  is  sufficiently 
small.  It  will  also  be  shown  that  the  instability  is  worse  if  too  many  lines  are  used. 

3.3  Analysis  of  Stability  of  MOL 


equation 


To  illustrate  the  stability  and  convergence  consider  the  problem  of  solving  Laplace’s 


=  0 


in  a  rectangular  domain  0  <  x  <  1,  -b  <  y  <  b,  with  the  following  Dirichlet  boundary  conditions: 


*<0,y)  -  *d,y)  =  0, 


i//(x,b)  =  ^(x,-b)  =  sin  7rx 


The  exact  solution  for  this  problem  is  known  to  be 


0(x,y)  = 


cosh  Try  sin  irx 
cosh  7rb 


We  now  consider  the  solution  of  this  problem  by  MOL.  Since  the  problem  contains  the  two 
lines  of  symmetry  x  =  xh  and  y  =  0,  we  can  reduce  the  region  of  interest  to  the  upper  left  quadrant  of 
the  rectangle,  i.e.  0  <  x  <  V2,  0  <  y  <  b.  N  -  1  interior  lines  are  drawn  parallel  to  the  x  axis  with  equal 
spacing  h  =  b/N,  so  that 


yn  =  nh  =  nb/N 


The  symmetry  conditions 


'  *  *  • _  _ 
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•M'^y)  =  o, 


(3.9) 


i//(x,-y)  =  i//(x,y), 


(3.10) 


are  applied. 

To  get  some  insight  into  the  stability  and  convergence  of  MOL  the  three-point  formula  (2.5) 
is  used  to  approximate  \pyy  in  (3.4)  giving 

-e  (^n+i  -  2^/n  +  i//n_t)/h2  =  0  (n  =  0,1,2,...,  N-  1)  (3.11) 


where  i//n(x)  is  the  approximation  for  i//(x,yn)  and  the  primes  indicate  differentiation  with  respect  to  x. 
To  the  system  (3.11)  we  add  the  appropriate  i.-jundary  and  symmetry  conditions. 


«^n(0)  =  0, 

(3.12) 

o' 

II 

3? 

.  e 
> 

(3.13) 

i//N  =  sin  jtx, 

(3.14) 

^-n(x)  =  *„(x),  n  =  1,2,...,  N-  1 

(3.15) 

3.3.1  Complementary  Function 

To  obtain  the  complete  solution  of  (3.11)-(3.15)  we  first  look  for  the  complementary 
function.  This  function  must  satisfy  (3.11),  (3.15)  and  also 


We  attempt  a  solution 


=  0 

iMx)  -  Cn  eMX 


and  substitute  into  (3.11)  to  give  the  recurrence  relation 

Cn+i  -  2z  Cn  +  Cn_j  =  0 

where 


(3.16) 


(3.17) 


•  -4sarv 
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For  the  complementary  function  to  satisfy  (3.15)  requires  C_j  =  C,  and  thus  using  (3.17)  we  obtain, 
for  C0  arbitrary, 


Ci  -  z  ■  Cq  -  Tj(z)C0 

C2  =  2zC,  -  C0  =  (2z2  -  1)C0  =  T2(z)C0 

C3  =  2zC2  -  Ci  =  2z(2z2  -  1)C0  -  zC0  =  (4z3  -  3z)C0  =  Tj(z)C0 


~  TN(z)C0 


(3.18) 


where  Tn(z)  is  the  Chebyshev  polynomial  of  order  n.  Now  to  satisfy  (3.16),  since  C0  f  0,  we  have 


0  =  Tn(z)  =  cosNO 


where  z  =  cosd.  Thus  8  takes  on  discrete  values  given  by  8m  where 


N0m  =  —  for  m  =  1,3, . .  . ,  2N  -  1 


,  ni  mir 

=  Zm  =  cos  8m  =  1  -  2  sin2  —  =  1  -  2  sin2  — 


and  so 


2  m7r 

"■J  ' 


(3.19) 


Taking  a  linear  combination  of  the  allowable  solution  and  using  the  fact  that  Tn(z)  -  cos  nd 
gives  us  the  complete  complementary  function 


*(ncV)  =  2  cos  ( Am  e"mX  +  Bm  e 

m=l,3 ...  ^ 


(3.20) 


for  constants  Am,  B„ 
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3.3.2  Particular  Integral 

In  this  case  we  follow  a  similar  procedure  to  that  of  finding  the  complementary  function 
except  that  we  attempt  a  solution 


=  Cn  sin  7rx 


(3.21) 


and  obtain  the  recurrence  relation 


C„+i  2z  Cn  +  Cn_|  -  0 


where  z  =  1  +  ir2h2/2.  Thus  we  obtain  the  same  relations  (3.18),  and  since  we  now  require  \J/N  -  sin  7rx 
we  have 


Cn  -  Tn(z)C0  -  1 


or 


C0 


1 

Tn<*) 


On  putting  cosh  0  =  z  and  using  the  fact  that  Tn(cosh  0)  =  cosh  n0,  the  particular  integral  (3.21) 
becomes 


<p)(x) 


cosh  n0 
cosh  N0 


sin  7rx 


(3.22) 


where  cosh  0  =  1  +  7r2h2/2. 

The  complete  solution  to  our  problem  is  then 

iMx)  =  «/'nP>(x)  +  ^nC)(x> 


3.3.3  Instability  in  the  Method  of  Lines 

Clearly  the  solution  we  require  by  MOL  is  the  particular  integral  (3.22).  The  complementary 
function  (3.20),  in  order  to  satisfy  the  boundary  conditions  (3.12)  and  (3.13),  requires  Am  =  Bm  =  0 
for  all  m.  However  in  solving  the  problem  numerically  it  may  be  that,  due  to  computer  round  off  and 
discretization  errors  in  the  numerical  integration,  the  Am  and  Bm  are  not  quite  zero.  In  this  case  we 

may  expect  a  large  growth  of  the  unwanted  solution  Am  e  01  in  (3.20)  since  /imx  may  grow  rather 
large.  The  size  of  the  unwanted  term  and  its  form  can  be  investigated  by  looking  at  the  largest  contri¬ 
bution  in  (3.20).  This  is  found  by  setting  m  =  2N  -  1  since  then 


. '!  - 
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2  .  /it  ir\  2  it 

um  =  —  sin - 1  =  —cos —  = 

m  h  \2  4N/  h  4N 


2N  t r 
—  cos  — 
b  4N 


Thus,  since  N  is  normally  greater  than  3,  nm  ~  2N/b  and  the  unwanted  growth  in  i£„c)(x)  is 


G2n_!  —  ^2N-l  cos  nit 


('-S> 


,2Nx/b 


The  cosine  factor  is  approximately  unity  in  magnitude  thus 


|G2N-l|  —  A2N-lexP^ 


2Nx\ 


b  / 


(3.23) 


It  can  be  seen  that  the  gradient  d/dx[i//„c,(x)]  will  grow  even  more  rapidly  since  this  will  be  given  by 


.  .  ,  2N  /2Nx\ 

|G2N-l|  -  a2N-i  b  exp  ^  b  J 


(3.24) 


In  our  particular  problem  we  really  require  ^J,(V&)  =  0  to  give  symmetry  about  x  =  V2  but  it  can  be  seen 
from  (3.24)  that  if  A2N_j  is  not  quite  zero  and  if  the  product  Nx  is  sufficiently  large  then  it  would  be 
impossible  to  satisfy  the  symmetry  condition. 

Now  due  to  the  truncation  errors  of  a  finite  word  length  machine  the  initial  conditions  at 
x  =  0  are  not,  in  general,  stored  exactly.  This  would  have  the  effect  of  Am  and  Bm  being  not  quite 
zero  and  so  lead  to  the  instability. 

Also,  even  with  an  infinite  word  length  machine,  the  discretization  errors  of  the  numerical 
integration  scheme  will  cause  instability.  This  can  be  seen  by  considering  the  first  step  from  x  =  0  to 
x  =  5x  using  the  exact  conditions  at  x  =  0  i.e.  ^n(0)  =  0  and  ^ „(0)  =  it  cosh  n0/cosh  N0.  We  obtain  a 
numerical  solution  at  x  =  Sx  given  by 


\pn  =  inexact)  +  en 
(exact)  +  5n 


where  en,  6n  are  the  errors  due  to  numerical  integration.  Thus  the  Am  and  Bm  now  satisfy 


2cos^(AmeMm6x  +  Bme'Mm6x)  =  en 


m 


2N 


2 

m 


niTlTT  UmSx 

Mrn  cos  —  (Am  e 


B  e  "m8x)  =  5 

Dm  c  f  wn 
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so  clearly  Am  and  Bm  are  not  zero  and  we  again  have  an  instability  problem. 

We  have  thus  shown  that  instability  is  unavoidable  but  we  do  not  know  how  serious  it  is 
until  we  attempt  solutions  to  practical  problems.  Since  we  now  know  the  form  of  the  unstable  term 
(3.23)  or  (3.24)  we  can  at  least  hope  to  minimize  its  effect. 

For  instance,  since  it  is  the  product  Nx  which  causes  instability,  it  may  in  some  cases  be 
worthwhile  using  smaller  integration  distances  x.  This  leads  to  the  idea  of  the  parallel  shooting 
technique  mentioned  at  the  end  of  Chapter  2.  Also  using  a  higher  (than  second)  order  difference 
scheme  for  the  i//yy  derivative  enables  us  to  use  fewer  lines  (N)  to  achieve  the  same  accuracy  for  this 
second  derivative.  In  turn  by  using  N  smaller  we  reduce  the  exponential  growth  in  (3.23)  or  (3.24). 

For  example  the  fourth  order  accurate  scheme  using  five  adjacent  lines  could  be  used.  This  scheme 
was  considered  in  Section  2.6  and  its  numerical  application  is  covered  in  Section  3.5.3. 

Notice  that  although  the  instability  is  unavoidable  it  is  no  worse  than  other  similar  techniques, 
used  particularly  in  the  Soviet  Union,  such  as  Telenin’s  method  (see  Gilinskii  et  al.  (1964))  or  the 
Method  of  Integral  Relations  (MIR).  Indeed  it  can  be  shown  that  the  instability  of  MIR  is  likely  to  be 
much  worse  than  that  of  MOL.  In  Appendix  A  the  MIR  instability  term  similar  to  (3.24)  is  derived 
and  shown  to  be 


A2N-i  exp 


(3.25) 


at  x  =  lh.  The  functions  (3.24)  and  (3.25),  without  the  A2n_j  ,  are  tabulated  in  Table  I  for  x  =  V2  and 
b  =  0.475  (the  value  used  in  our  example  later).  From  this  table  it  can  be  seen  that  instability  of  MIR 
is  significantly  worse. 

Even  so  MIR  has  been  used  successfully  in  the  Soviet  Union  and  also  by  Holt  (1977)  for 
many  years.  In  Chapter  5  of  Holt’s  book  application  of  MIR  to  several  problems  such  as  the  supersonic 
blunt  body  problem  and  the  laminar  boundary  layer  equations  is  explained. 


3.4  Convergence  of  the  Particular  Integral 

It  is  of  interest  to  compare  the  particular  integral  (3.22)  with  the  exact  solution  (3.7).  This 
is  accomplished  by  expanding  the  inverse  hyperbolic  function  in  0  =  cosh-1  (1  +  ir2h2/2)  giving 


6  =  jrh  jl  -  ^j-+  0(h4) 


and  then  expanding  the  particular  integral  in  a  Taylor  series  about  h  =  0 


*n(x)  - 


cosh  iryn 
cosh  7rb 


■  sin  irx 


ir3b3  (  --  ,  4  . 

1+ -  Itanh  7rb - tanh  wynJ  +  0(N  4)| 

24N2  '  b 


yn 


f„)  +  0(N-4)j 


Thus  the  particular  integral  converges  to  the  exact  solution  (3.7)  and  has  error  0(N  2).  This  error  is 

expected  since  we  use  a  formula  of  0(h2)  accuracy  for  the  second  derivative  ^yv  and  h  =  77. 

3  N 
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3.5  A  Linear  Example  to  Illustrate  the  Stability  and  Convergence 

In  this  section  we  solve  numerically  the  problem  posed  in  Section  3.3  and  given  by  Equations 
(3.4)-(3.6)  using  b  =  0.475.  The  exact  solution  (3.7)  is  tabulated  in  Table  II  for  values  of  x  =  0.125, 
0.25,  0.375,  0.5  and  for  values  of  y  =  0,  0.475/3  and  2/3  X  0.475;  also  shown  are  values  of  i//x(x  =  0) 
at  the  above  y  values. 

We  will  use  the  three  line  difference  scheme  (2.5)  in  order  to  get  a  better  understanding  of 
the  stability  since  in  this  case  we  know  the  analytic  MOL  solution  is  given  by  (3.22).  With  the  three 
line  scheme  we  will  investigate  the  effects  of  parallel  shooting,  of  varying  the  number  of  lines,  and  of 
varying  the  integration  step  size  6x. 

We  will  then  use  the  five  line  scheme  and  show  that  accurate  solutions  can  be  obtained  using 
only  a  few  lines. 

3.5.1  The  Three  Line  Analytical  Solution 

Let  us  consider  the  analytical  solution  to  the  MOL  equations  which  is  given  by  the  particular 
integral  (3.22)  of  the  ordinary  differential  Equations  (3.11).  This  analytical  solution  represents  the 
best  we  can  hope  to  achieve  numerically  since  it  does  not  contain  any  round  off  or  discretization  errors 
and  hence  the  complementary  function  (3.20)  is  zero.  Remember  that,  numerically,  we  expect  (3.20) 
to  be  non  zero  and  hence  to  cause  some  inaccuracy.  The  extent  of  the  inaccuracy  is  what  we  want  to 
observe.  Thus  comparison  of  our  numerical  solution  with  the  analytical  MOL  solution,  rather  than 
with  the  exact  solution  (3.7),  is  more  meaningful  at  this  stage.  The  analytical  MOL  solution  is  given  in 
Table  III  at  the  same  x  and  y  values  as  listed  for  Table  II;  we  present  results  for  N  =  3,6,9,12  and  15 
(i.e.  the  number  of  lines  used  to  divide  the  region  0  <  y  <  b  into  strips). 

3.5.2  The  Three  Line  Numerical  Scheme 

This  section  is  split  into  three  parts  which  illustrate  three  different  ways  of  overcoming  the 
instability. 

The  present  results  are  obtained  using  the  Runge  Kutta  Gill  integration  method  which  has 
truncation  error  0(6x5 ).  One  step  of  the  Newton  iteration  scheme  (2.10)  is  used  to  find  the  missing 
initial  conditions.  In  the  tables  to  follow  we  print  values  of  \p  and  the  initial  slopes  i^x(x  =  0)  at  the 
(x,y)  values  mentioned  in  Section  (3.5).  We  also  print  N  (the  number  of  lines),  the  integration  step 
size  DELX,  and  the  matching  point  XSHOOT  used  in  the  parallel  shooting  technique.  If  XSHOOT  =  0.5 
then  only  one  shooting  range  is  used  and  we  print  out  the  SUM  SQUARES  GRADIENTS  which  repre¬ 
sents  2p^  at  x  =  0.5  where  p,  =  9^ /9x;  because  of  symmetry  about  x  =  0.5  this  quantity  should  be 
zero.  If  XSHOOT  f  0.5  then  two  shooting  ranges  have  been  used  as  shown  below. 

(a)  Effect  of  Varying  the  Length  of  Shooting  (Parallel  Shooting) 

In  order  to  help  to  overcome  the  instability  we  can  use  a  parallel  shooting  technique  as 

shown  diagrammatically  below. 


u 

0 


x  (shoot) 


x 


We  first  estimate  Pj(0),  integrate  the  differential  equations  from  0  to  x(shoot)  and  obtain 
values  \pf  and  p?,  say,  on  line  i  at  x  =  x(shoot).  Also  using  estimates  of  and  ps  at 
x  =  x( shoot)  (designated  and  p?H  say)  we  integrate  the  ordinary  differential  equations 
to  the  end  of  the  range  x  =  0.5  where  we  obtain  values  of  pj  =  Pj(0.5). 


f 

t 


-Avar-. 
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Be  suitably  choosing  Pj(0),  and  p?H  we  hope  to  solve  the  system  of  equations 


Pj  (0.5)  =  0  (symmetry  about  x  =  0.5) 

-  *?«  =  0 

matching  at  x  =  x(shoot) 

pf  -  pfH  =  0  (3.26) 


The  optimum  choice  of  Pj(0),  and  p?H  is  made  using  Newton’s  method  (2.10).  Our 
linear  illustrative  problem  remains  linear  since  the  matching  conditions  at  x  =  x( shoot)  are 
linear.  Thus  one  step  in  the  above  method  is  all  that  is  required. 

Table  IV  shows  the  results  of  varying  x( shoot)  with  N  =  15  and  fix  =  1/32.  It  can  be  seen 
that  instability  is  least  when  using  the  smallest  shooting  range  (i.e.  x( shoot)  =  0.25).  Shoot¬ 
ing  the  whole  range  is  clearly  hopeless  while  using  x( shoot)  =  0.125  or  0.375  gives  fairly 
reasonable  results.  Using  x( shoot)  =  0.25  gives  good  stability  and  hence  a  solution  accurate 
to  almost  4  figures  compared  to  the  analytic  MOL  solution  given  in  Table  III.  The  quantity 
SUM  SQUARES  RESIDUALS  printed  in  the  tables  represents  the  sum  of  squares  of  the  left 
hand  sides  of  (3.26).  The  two  entries  given  for  x  =  x(shoot)  list  values  of  i//?  (obtained  after 
integration  from  zero)  and  i//?H  (the  best  values  of  ^(XSHOOT)  to  minimize  the  RESIDUALS 
of  (3.26)). 

(b)  Effect  of  Varying  the  Integration  Step  Length  fix 

It  has  been  shown  that  the  instability  is  caused  by  terms  of  the  type  exp(p2N-j  *)•  °n 
integrating  numerically  by  the  Runge  Kutta  Gill  method  we  find  that  after  k  steps  the  term 
exp(p2N-i  k  Sx)  is  represented  by  Ek  where 


A2  A3 
E  =  1  +  A  +  —  +  — 
2!  3! 


+ 


A^ 

4! 


and  A  =  p2w_1  Sx.  Thus  the  error  grows  in  the  numerical  integration  like  Ek  rather  than 
exp(kA).  Ek  is  always  less  than  exp(kA)  for  A  >  0  and  the  ratio  Ek/exp(kA)  shrinks  rapidly 
with  increasing  A.  For  example  when  fix  =  1/4  (i.e.  k  =  2  for  integration  over  the  whole 
range)  Ek/exp(kA)  is  about  10-4  for  N  =  12. 

Thus  we  may  expect  instability  to  improve  as  Ax  is  increased  (but  of  course  the  particular 
integral  would  not  be  so  accurately  represented). 

Table  V  shows  the  improvement  in  stability  for  N  =  15  and  x( shoot)  =  0.5.  We  can  see  that 
we  progress  from  the  unstable  cases  of  fix  =  1/32  or  1/16  to  a  fairly  stable  solution  with 
fix  =  1/8  and  to  a  stable  solution  with  fix  =  1/4.  Referring  to  accuracy  the  latter  two 
solutions  are  reasonably  good  except  at  the  point  x  =  1/2,  y  =  2/3*  0.475. 

(c)  Effect  of  Varying  the  Number  of  Lines 

It  can  be  seen  from  Table  VI  that  the  effect  of  decreasing  the  number  of  lines  improves 
enormously  the  behaviour  of  the  instability.  The  N  =  15  solution  is  hopelessly  unstable, 
while  using  N  =  12  gives  a  mildly  unstable  situation  (the  results  are  not  accurate  compared 
to  the  analytical  MOL  solution  of  Table  III).  Using  N  =  9,  6  or  3  stability  does  not  seem  to 


-27  - 


be  any  problem  and  the  accuracy  of  the  solution  is  very  good  compared  to  the  analytic  MOL 
solution  (although  not  necessarily  compared  to  the  exact  solution  given  in  Table  II). 

To  get  good  accuracy  compared  to  the  exact  solution  using  the  3  line  difference  scheme  we 
would  clearly  have  to  use  many  lines  (N  =  57  would  give  4  decimal  figures  accuracy)  but 
this,  numerically,  would  give  instability.  Thus  we  clearly  must  use  a  higher  order  difference 
scheme  to  represent  the  derivatives  so  that  we  can  minimize  the  number  of  lines  necessary  to 
achieve  an  accurate  representation  of  the  particular  integral.  This  scheme  is  now  investigated 
numerically. 

3.5.3  The  Five  Line  Numerical  Scheme 

As  we  saw  from  the  three  line  scheme  applied  to  our  linear  example  we  can  expect  good 
stability  with  N  <9  but  the  accuracy  may  not  be  sufficiently  accurate.  To  keep  the  stability  and  also 
get  accuracy  we  consider  the  fourth  order  accurate  scheme  using  (2.17b)  to  represent  the  derivative 
\pyy  (with  (2.20)  used  next  to  the  boundary  line  y  =  b).  This  time  we  only  consider  shooting  the  one 
range  0  <  x  <  0.5  since  sufficient  accuracy  and  stability  is  obtained  in  this  case  (using  N  <  9), 

Table  VII  lists  values  of  i//  and  the  initial  slopes  i|/x(x  =  0)  for  some  values  of  N  and  6x.  It 
can  be  seen,  by  comparing  with  Table  II,  that  almost  4  significant  figure  accuracy  is  obtained  using 
N  =  6  and  Sx  =  0.0625.  Notice  also  that  using  the  coarsest  scheme,  N  =  3  and  6x  =  0.25,  produces 
results  with  at  worst  2%  error  in  \p  and  1%  error  in  the  initial  slopes.  The  N  =  12  results  again  indicate 
instability  when  Sx  is  too  small. 

In  terms  of  computer  time  the  method  is  very  efficient  as  can  be  seen  from  the  N  =  6, 

Sx  =  0.0625  solution.  This  required  only  0.3  seconds  on  an  IBM  3032  computer.  The  efficiency  is 
confirmed  by  further  computer  times  quoted  in  the  nonlinear  examples  covered  in  the  next  chapter. 

3.6  Summary 

The  three  line  scheme  (2.5)  to  represent  the  derivative  has  been  used  mainly  to  illustrate  the 
instability  inherent  in  the  method  of  lines.  It  has  been  confirmed  that  three  techniques  help  to  over¬ 
come  the  instability  but  the  latter  two  decrease  the  accuracy  of  the  particular  integral  which  we  seek 
to  obtain.  It  has  been  shown  that,  for  our  linear  example,  the  higher  order  difference  scheme  using 
five  adjacent  lines  (2.17)  gives  an  accurate  representation  of  the  particular  integral  while  at  the  same 
time  giving  a  stable  solution.  Thus  use  of  the  five  line  scheme  is  always  recommended. 

Notice  also  that  three  line  schemes  of  order  (N-4)  accuracy  in  the  particular  integral  may  be 
used  in  certain  cases.  These  schemes  are  given  in  Appendix  B;  they  require  calculating  the  derivatives 
Pi  in  the  ordinary  differential  equations  by  solution  of  a  tridiagonal  system  of  equations.  This  can  be 
achieved  efficiently  by  the  Thomas  algorithm  (von  Rosenburg,  1969). 
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TABLE I 

TABLE  OF  THE  FUNCTIONS  (3.24)  AND  (3.25)  USING  b  =  0.475 
TO  ILLUSTRATE  THE  INSTABILITY  OF  MOL  AND  OF  THE 
METHOD  OF  INTEGRAL  RELATIONS  FOR  N  TOO  LARGE 


N 

(3.24) 

(3.25) 

2 

5.7,2 

4.5,4  (=  4.5  X  104) 

3 

7.0,3 

3.0, 10 

4 

7.6,  4 

4.2, 18 

5 

7.8,  5 

1.3,29 

6 

7.7,6 

8.1,41 

9 

6.4,9 

>107S 

12 

4.7, 12 

>1075 

TABLE  II 

EXACT  SOLUTION  TO  MOL  LINEAR  PROBLEM 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1638 

0.1845 

0.2518 

0.250 

0.3027 

0.3409 

0.4653 

0.375 

0.3955 

0.4454 

0.6079 

0.500 

0.4281 

0.4821 

0.6580 

INITIAL  SLOPES 

1.3449 

1.5147 

2.0671 

J 


TABLE  III 


ANALYTICAL  MOL  SOLUTIONS  GIVEN  BY  3.22 


N  =  3 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1660 

0.1866 

0.2533 

0.250 

0.3068 

0.3448 

0.4680 

0.375 

0.4009 

0.4505 

0.6115 

0.500 

0.4339 

0.4876 

0.6619 

INITIAL  SLOPES 

1.3632 

1.5318 

2.0794 

N  =  6 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1644 

0.1850 

0.2522 

0.250 

0.3037 

0.3419 

0.4660 

0.375 

0.3969 

0.4467 

0.6088 

0.500 

0.4296 

0.4835 

0.6590 

INITIAL  SLOPES 

1.3495 

1.5190 

2.0702 

N  =  9 


X 

Y  ■  0 

0.475/3 

2*0.475/3 

0.125 

0.1641 

0.1847 

0.2520 

0.250 

0.3032 

0.3414 

0.4656 

0.375 

0.3961 

0.4460 

0.6083 

0.500 

0.4287 

0.4828 

0.6584 

INITIAL  SLOPES 

1.3469 

1.5166 

2.0685 

N  =  12 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1640 

0.1846 

0.2519 

0.250 

0.3030 

0.3412 

0.4654 

0.375 

0.3958 

0.4458 

0.6081 

0.500 

0.4285 

0.4825 

0.6582 

INITIAL  SLOPES 

1.3460 

1.5158 

2.0679 

N  =  15 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1639 

0.1846 

0.2519 

0.250 

0.3029 

0.3411 

0.4654 

0.375 

0.3957 

0.4456 

0.6080 

0.500 

0.4283 

0.4824 

0.6581 

INITIAL  SLOPES 

1.3456 

1.5154 

2.0676 

TABLE  IV 


EFFECT  OF  VARYING  X(SHOOT) 


N  =  15  DELX  =  0.03125  XSHOOT  =  0.125 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1635 

0.1840 

0.2511 

0.125 

0.1635 

0.1840 

0.2511 

0.250 

0.3023 

0.3401 

0.4630 

0.375 

0.3955 

0.4452 

0.6015 

0.500 

0.4284 

0.4833 

0.6315 

INITIAL  SLOPES 

1.3422 

1.5108 

2.0623 

SUM  SQUARES  RESIDUALS  6.885D  02 


N  =  25  DELX  =  0.03125  XSHOOT  =  0.250 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1639 

0.1846 

0.2519 

0.250 

0.3029 

0.3411 

0.4654 

0.250 

0.3029 

0.3411 

0.4654 

0.375 

0.3958 

0.4457 

0.6081 

0.500 

0.4284 

0.4824 

0.6582 

INITIAL  SLOPES 

1.3457 

1.5155 

2.0677 

SUM  SQUARES  RESIDUALS  1.609D-04 


♦  N  =  25  DELX  =  0.03125  XSHOOT  =  0.375 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1657 

0.1860 

0.2526 

0.250 

0.3072 

0.3441 

0.4667 

0.375 

0.4028 

0.4701 

0.6643 

0.375 

0.4075 

0.4500 

0.6099 

0.500 

0.4354 

0.4870 

0.6601 

INITIAL  SLOPES 

1.3595 

1.5265 

2.0732 

SUM  SQUARES  RESIDUALS  6.799D  01 


N  =  15  DELX  =  0.03125  XSHOOT  =  0.500 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.0363 

0.0805 

0.1980 

0.250 

0.0129 

0.1170 

0.3574 

0.375 

-0.1656 

0.0768 

0.4544 

0.500 

0.6212 

-8.3976 

-1.2045 

INITIAL  SLOPES 

0.3648 

0.7049 

1.6394 

SUM  SQUARES  GRADIENTS  1.050D  06 


TABLE  V 


EFFECT  OF  VARYING  6x 


N  =  15  DELX  =  0.06250  XSHOOT  =  0.500 


X 

Y  =  0 

WHEEBZMm 

2*0.475/3 

0.125 

0.1583 

0.1812 

0.2510 

0.250 

0.2873 

0.3347 

0.4649 

0.375 

0.3473 

0.4450 

0.6011 

0.500 

0.5551 

1.0065 

0.0971 

INITIAL  SLOPES 

1.3044 

1.4881 

2.0591 

SUM  SQUARES  GRADIENTS  3.905D  03 


N  =  15  DELX  =  0.12500  XSHOOT  =  0.500 


X 

Y  *  0 

0.475/3 

2*0.475/3 

0.125 

0.1636 

0.1843 

0.2516 

0.250 

0.3022 

0.3406 

0.4652 

0.375 

0.3944 

0.4448 

0.6249 

0.500 

0.4260 

0.4838 

7.1353 

INITIAL  SLOPES 

1.3435 

1.5136 

2.0653 

SUM  SQUARES  GRADIENTS  8.146D-02 


N  =  15  DELX  =  0.25000  XSHOOT  =  0.500 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.250 

0.3035 

0.3413 

0.4638 

0.500 

0.4302 

0.4874 

0.7246 

INITIAL  SLOPES 

1.3500 

1.5165 

2.0533 

SUM  SQUARES  GRADIENTS  7.961D-07 


TABLE  VI 


EFFECT  OF  VARYING  THE  NUMBER  OF  LINES 


N  =  15  DELX  =  0.03125  XSHOOT  =  0.500 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.0363 

0.0805 

0.1980 

0.250 

0.0129 

0.1170 

0.3574 

0.375 

-0.1656 

0.0768 

0.4544 

0.500 

0.6212 

-8.3976 

-1.2045 

INITIAL  SLOPES 

0.3648 

0.7049 

1.6394 

SUM  SQUARES  GRADIENTS  1.050D  06 


N  =  12  '  DELX  =  0.03125  XSHOOT  =  0.500 


X 

o 

II 

> 

0.475/3 

2*0.475/3 

0.125 

0.1612 

0.1824 

0.2507 

0.250 

0.2967 

0.3363 

0.4631 

0.375 

0.3836 

0.4381 

0.6047 

0.500 

0.3944 

0.5251 

0.6659 

INITIAL  SLOPES 

1.3249 

1.4983 

2.0586 

SUM  SQUARES  GRADIENTS  2.130D  01 


N  =  9  DELX  =  0.03125  XSHOOT  =  0.500 


X 

Y  =■  0 

0.475/3 

2*0.475/3 

0.125 

0.1641 

0.1848 

0.2520 

0.250 

0.3032 

0.3414 

0.4656 

0.375 

0.3962 

0.4460 

0.6083 

0.500 

0.4289 

0.4826 

0.6585 

INITIAL  SLOPES 
_ 

1.3470 

1.5167 

2.0685 

SUM  SQUARES  GRADIENTS  2.326D-04 


N  =  6  DELX  =  0.03125  XSHOOT  =  0.500 


X 

Y  =  0 

0.475/3 

2*0.475/3 

0.125 

0.1644 

0.1850 

0.2522 

0.250 

0.3037 

0.3419 

0.4660 

0.375 

0.3969 

0.4467 

0.6088 

0.500 

0.4296 

0.4835 

0.6590 

INITIAL  SLOPES 

1.3495 

1.5190 

2.0702 

SUM  SQUARES  GRADIENTS  3.526D-09 


N  =  3  DELX  =  0.03125  XSHOOT  -  0.500 


X 

-< 

II 

o 

0.475/3 

2*0.475/3 

0.125 

0.1660 

0.1866 

0.2533 

0.250 

0.3068 

0.3448 

0.4680 

0.375 

0.4009 

0.4505 

0.6115 

0.500 

0.4339 

0.4876 

0.6619 

INITIAL  SLOPES 

1.3632 

1.5318 

2.0794 

SUM  SQUARES  GRADIENTS  1.841D-16 


NUMERICAL  SOLUTION  OF  LINEAR  PROBLEM  BY  THE  FIVE  LINE  SCHEME 
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CHAPTER  4.0  TWO  NONLINEAR  EXAMPLES 

4.1  Introduction 

The  linear  problem  covered  in  the  last  chapter  was  relatively  easy  to  solve  since  only  one 
step  in  the  Newton  scheme  was  required  tc  obtain  a  solution.  In  nonlinear  cases  we  have  to  be  much 
more  careful  about  initial  estimates  and  we  must  be  sure  of  using  an  efficient  numerical  optimisation 
scheme  in  order  to  find  the  optimum  unknown  values.  It  is  believed  that  the  method  due  to  Powell 
(1965)  is  one  of  the  most  efficient  methods  (see  Section  2.3  for  a  further  coverage  of  Powell’s  method) 
for  minimizing  a  sum  of  squares  of  nonlinear  functions  and  as  such  is  used  in  most  of  our  computations 
presented  here,  the  exception  being  the  delta  wing  computations  of  Klunker  et  al.  (1971)  who  used 
the  modified  Newton  scheme  mentioned  in  Section  2.2.2. 

The  first  example  is  a  structural  problem  and  as  such  has  been  solved  by  the  finite  element 
method  (Dixon  1971);  MOL  results  are  compared  with  the  finite  element  solution. 

The  second  example  is  a  conical  flow  problem  encountered  in  aerodynamics.  We  solve  the 
full  Euler  equations  comprising  five  simultaneous  partial  differential  equations  and  also  determine  the 
location  of  the  shock  wave  attached  to  the  cone  apex  (see  Fig.  4.3).  A  finite  difference  solution  of  the 
elliptic  equations  would  be  very  time  consuming  and  the  earlier  methods  of  solving  this  problem  con¬ 
sidered  the  complete  hyperbolic  system  using  the  cone  axis  as  the  time-like  direction.  Comparisons  of 
the  MOL  solution  to  the  hyperbolic  marching  method  (Babenko  1966;  Gonidou  1968)  are  made.  It  is 
found  that  MOL  solutions  are  some  fifty  times  quicker  while  obtaining  results  of  equal  accuracy.  The 
MOL  results  for  a  circular  cone  have  been  tabulated  in  AGARDograph  137  (Jones  1969)  and  have 
served  as  standard  reference  tables  for  inviscid  flow  about  circular  cones. 

Finally  some  results  are  shown  from  MOL  solutions  for  flow  about  the  compression  side  of 
a  delta  wing  at  incidence  to  a  supersonic  stream.  Again  a  comparison  is  made  with  a  hyperbolic  march¬ 
ing  method. 

Other  nonlinear  examples  are  discussed  in  Jones,  South  and  Klunker  (1972). 

4.2  The  Simply  Supported  Square  Plate 

For  this  example  we  use  the  following  notation 

constant  =  tE/(l-i>2) 
flexural  rigidity,  Et3/12  (l-»>2) 
modulus  of  elasticity,  =  2  X  1010 
plate  thickness,  =  0.002 
number  of  strips  used,  h  =  '/2/N 

=  ux 

=  vx 

normal  pressure.  Varies  from  0  to  80 

=  wx 

“  01 X 

mid  plane  displacements  in  x  and  y  directions 
transverse  displacement 
substituted  for  wxx  +  wyy 

step  length  used  for  integrating  the  MOL  ordinary  differential  equations 
strip  width  in  y  direction 
Poisson’s  ratio,  =  0.3 
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The  example  is  a  structural  problem  which  is  defined  mathematically  by  a  set  of  rather 
complicated  nonlinear  partial  differential  equations.  These  are 

ux x  +  wxwxx  +  v(vxy  +  wy  wxy)  + 

*(uyy  +  vxy  +  wx  wyy  +  wy  wxy)  =  0  (4.1a) 

Vyy  +  Wy  Wyy  +  H  Ujj  y  +  Wj^W^)  +  'A(l  ~  V) 

•(vxx  +  uxy  +  Wy  wxx  +  wx  wxy)  =  0  (4.1b) 


tE 

DV4W  =  Q  +  -  |[UX  +  ‘/2WX  +  l>(Vy  +  '/2Wy)] 

1-v 2 

•Wxx  +  [Vy  +  '/2Wy  +  V(UX  +  ‘/2WX)]  Wyy 

+  (l-V)(Uy  +  VX  +  WxWy)Wxy|  (4-lC) 

subject  to  boundary  conditions  on  the  simply  supported  square  plate  o<x<l,o<y<l 

x  =  o:  u  =  v  =  w  =  wxx  =  o  (4.2a) 

y  =  o:  u  =  v  =  w  =  wyy  =  o  (4.2b) 


with  lines  of  symmetry  x  =  V2,  y  =  Vi. 


To  use  MOL,  Equations  (4.1)  are  written  as  differential  equations  of  first  order  in  x  which  is 
the  variable  to  be  left  continuous.  We  also  for  convenience,  introduce  a  variable  a  =  V2w  (hence,  on 
the  boundaries,  a  =  o),  and  suitably  scale  the  dependent  variables  u  =  Cu,  v  =  Cv,  w  =  C,/jw,  a  =  C1/ja. 
Then,  dropping  the  bar  notation,  the  Equations  (4.1)  can  be  written  as  the  system  of  equations 


ux  =  p  (4.3a) 

vx  =  q  (4.3b) 


w 


X 


r 


(4.3c) 


r 
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ax  =  s  (4.3d) 


=  a  -  w 


yy 


(4.3e) 


px  =  -  [rrx  +  v( qy  +  wy  ry)  +  Vi(l-v) 

"(Uyy  +  qy  +  rwyy  +  wyry)]  (4.3f) 

qx  =  -  2(l-v)~'  [vyy  +  wy  wyy  +  i>(Py  +  r ry)]  -  py  -  wy  rx  -  rry  (4.3g) 

Dsx  =  -  Dayy  +  C'/2Q  +  [p  +  'Ax2  +  v(\y  +  V&w2)]rx 

+  [vy  +  '/2  Wy  +  i»(p  +  l^r2)]Wyy  +  (l^^Uy  +  Q  +  FW^Ty  (4.3h) 


Now  Equations  (4.3)  are  to  be  satisfied  on  each  dividing  line  (except  at  y  =  o)  which  are 
numbered  0,1,2  ...  N  say  corresponding  to  increasing  y.  Thus  y0  =  0,  yN  =  ‘/z  and  the  strip  width 
h  =  1/2N.  Symmetry  on  y  =  lh  can  be  conserved  by  introducing  image  lines  at  lh  +  h  and  Vs  +  2h  and 
these  lines  are  numbered  N  +  1  and  N  +  2.  Then  for  symmetry 


vN+k - vN-k 


^N+k  “  ^N-k 


for  k  =  1,2  while  for  the  other  dependent  variables  the  equality  holds  (e.g.  wN+k  =  wN_k). 

The  partial  derivatives  with  respect  to  y  occurring  in  Equations  (4.3)  are  replaced  by  the  five 
line  difference  formulas  (2.17)  for  i  =  2,3  . . .  N,  while  on  line  1  we  use  formulas  (2.19)  and  (2.20). 

Now  we  can  apply  MOL  by  estimating  p,q,r  and  s  at  each  line  i  at  x  =  o.  These  estimates  are 
made  by  first  obtaining  a  solution  with  Q  near  zero  since  at  Q  =  0  we  know  p  =  q  =  r  =  s  =  o.  We  then 
increase  Q  gradually  and  obtain  initial  estimates  of  p,q,r  and  s  at  x  =  o  by  extrapolation  from  previous 
solutions;  in  our  case  quadratic  extrapolation  was  used  as  soon  as  possible.  We  then  improve  the 
estimates  using  Powell’s  method  in  order  that 


r2  +  q2  +  u2  +  s2 


(4.4) 


at  each  line  at  x  =  xh  is  minimized.  The  latter  must  be  minimized  since  we  require 


r  =  q  =  u  =  s  =  o 
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at  x  =  xh  to  preserve  symmetry  there.  In  the  practical  computation  it  turned  out  that  weighting  factors 
had  to  be  applied  in  (4.4)  so  that  each  term  was  roughly  of  the  same  order. 

The  integration  method  used  was  Hamming’s  (1959)  predictor  modifier  corrector  scheme  of 
0(5x5)  accuracy.  Starting  with  N  =  3  various  step  lengths  5x  were  used  i.e.  6x  =  'A/m  where  m  =  3,4,5, 
6,7,8,10.  This  procedure  indicated  the  step  length  8x  needed  to  give  sufficient  accuracy.  It  was  found 
that  5x  =  1/20  was  sufficiently  small  since  these  results  differed  only  a  little  from  6x  =  1/16. 

We  then  used  N  =  4  lines  (5x  =  1/20)  and  compared  the  results  with  N  =  3  for  Q  =  32.  Some 
quantities  were  not  sufficiently  accurate,  e.g.  wyy,  so  we  next  used  N  =  5  and  again  compared  results  at 
Q  =  32.  Some  of  these  comparisons  are  shown  in  Figures  4.1.  It  was  found  that,  with  N  =  5,  solutions 
were  rather  difficult  to  obtain  for  Q  >  32;  this  was  due  to  the  instability  becoming  more  significant  as 
the  nonlinearity  of  the  problem  increased  (parallel  shooting  (Section  2.9)  would  have  to  be  used  to 
obtain  solutions  with  higher  Q).  However  using  the  one  shooting  range  we  obtained  satisfactory  results 
with  N  =  3  or  4.  Solutions  at  Q  =  80  are  compared  in  Figures  4.2  with  a  9  element  finite  element 
solution  of  Dixon  (1971).  It  is  seen  that  good  agreement  is  obtained. 

The  time  taken  for  the  MOL  solutions  for  some  40  values  of  Q  between  0  and  80  was  as 

follows 


N  =  3, 


2.5  mins. 


N  =  4, 


3.5  mins. 


on  an  IBM  3032  computer. 

4.3  Conical  Flow  Problems 

The  method  of  lines  seems  to  be  ideally  suited  to  problems  of  this  type  since  it  is  not  clear 
how  one  would  solve  the  elliptic  system  by  finite  difference  methods.  The  previous  finite  difference 
methods  (e.g.  Babenko  1966)  solved  the  hyperbolic  system  of  equations  by  a  marching  technique  using 
the  cone  axis  z  (see  Fig.  4.3)  as  the  time-like  direction.  But  these  marching  methods  were  rather  slow, 
for  example,  Gonidou  (1969)  reported  one  hour  computation  for  a  certain  configuration  and,  as  many 
solutions  for  different  configurations  were  required,  it  seemed  that  a  more  efficient  computational 
technique  was  needed.  MOL  turned  out  to  be  very  efficient  with  between  5  and  20  seconds  needed  for 
a  solution.  For  flows  about  circular  cones  a  set  of  tables  of  MOL  results  has  been  published  by  Jones 
(1969);  results  for  some  1200  configurations  are  listed  there. 

The  physical  problem  is  shown  in  Figure  4.3.  In  supersonic  flow  there  is  a  shock  wave 
attached  to  the  tip  of  the  conical  body  and  the  flow  field  behind  the  shock  is  such  that  quantities 
along  any  ray  emanating  from  the  tip  are  constant.  We  wish  to  find  the  shock  wave  shape  and  the  flow 
field  variables  between  the  conical  body  in  supersonic  flow  and  its  attached  shock  wave.  A  cylindrical 
co-ordinate  system  (z,r,0)  is  adopted  with  the  z  axis  along  the  axis  of  the  conical  body  (which  may  for 
example  be  a  circular  cone). 

The  equation  of  the  given  body  can  be  written  in  the  form,  r  =  zG(0)  say,  and  we  let  the 
equation  of  the  unknown  attached  shock  wave  be  r  =  zF(0).  The  full  three-dimensional  equations  of 
motion  (momentum,  continuity  and  energy  conservation)  can  be  written  in  matrix  form  as 


,  9X 
dz 


o-  dX 

+  B  —  + 
9r 


+ 


D' 


=  0, 


(4.5) 
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where  A'B'C'  are  (5  X  5)  matrices,  D'  is  a  column  vector  and  X  is  also  a  column  vector  given  by 


X  = 


where  u,v,w  are  the  velocity  components  in  the  (z,r ,0)  directions,  respectively,  p  is  the  pressure  and  p 
the  density.  The  matrices  and  the  vector  D'  consist  of  elements  which  are  functions  of  u,v,w,p,  and  p; 
their  exact  form  can  be  found  in  Jones  (1968).  A  cross  section  (z  =  const)  of  the  flow  field  is  shown 
in  Figure  4.4a;  here  the  flow  and  body  are  assumed  to  be  symmetrical  about  0  =  0,  ir  as  is  usually  the 
case.  The  boundary  conditions  to  be  satisfied  are  the  Rankine-Hugoniot  relations  at  the  shock  wave 
which  can  be  written  in  the  form 


X  =  f(a,7>Ioo,0,F,F'), 


(4.6) 


where  f  is  a  column  vector  whose  elements  are  functions  of  the  listed  arguments.  The  first  three 
arguments  in  (4.6)  are  a  the  angle  of  incidence  which  is  the  angle  that  the  directior  of  the  free  stream 
makes  with  the  z  axis,  7  the  ratio  of  specific  heats  and  the  free-stream  Mach  number.  Since  these 
three  arguments  are  known  for  a  given  problem  it  follows  that  the  elements  of  X  are  known  at  the 
shock  once  the  equation  of  the  shock  r  =  zF(0)  is  known.  F'  is  found  from  (2.17a). 

On  the  body  the  normal  velocity  should  be  zero  and  can  be  written 


uG  -  v  +  (l/G)(dG/d0)w  =  0.  (4.7) 

Now  it  is  known  that  the  equations  of  motion  (4.5)  can  be  reduced  to  two  dimensions  since 
the  flow  is  conical.  A  suitable  transformation  to  do  this  and  also  one  which  makes  the  boundaries 
easier  to  handle  is  given  by 


x  =  z, 

t  -  [r  -  zG(0)]/z[F(0)  -  G(0)), 

0  =  0. 


It  is  now  seen  that  the  body  r  =  zG(0)  and  shock  wave  r  =  zF(0)  are  transformed  to  the  lines  £  =  0  and 
£  =  1,  respectively,  see  Figure  4.4b.  The  equations  of  motion  (4.5)  are  transformed  to 


B(3X/6£)  +  C(9  X/30)  +  D  =  0, 


(4.8) 


where  B,C  are  (5  X  5)  matrices  and  D  is  a  column  vector.  The  term  dX/dx  is  omitted  from  the  above 
equation  since  the  flow  is  conical  and  3X/3x  is  zero.  Now  we  can  consider  the  equations  at  unit 
distance  x  *  z  =  1.  Hence  the  problem  is  reduced  to  that  of  finding  solutions  of  (4.8)  in  the  region 
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O<f<l,O<0<;r  (assuming  symmetry)  subject  to  boundary  conditions  (4.6)  at  £  =  1  and  (4.7)  at 

S  =  o. 


The  method  of  lines  is  applied  in  the  (£,0)  plane;  symmetry  conditions  at  0  =  0  and  at  0  =  jt 
are  satisfied  by  introducing  image  lines  in  the  usual  manner.  An  estimate  for  F(0)  at  each  of  the  lines 
is  made,  F'(0)  is  obtained  from  (2.17a)  and  substitution  in  (4.6)  gives  X(£  =  1).  Equations  (4.8)  are 
next  reduced  to  ordinary  differential  equations  by  writing  3X/S0  at  each  line  in  the  finite-difference 
form  (2.17a).  Integration  of  these  ordinary  differential  equations  is  then  made  from  {  =  1  to  £  =  0 
where  Equation  (4.7)  must  be  satisfied  at  each  dividing  line.  The  shock  shape,  i.e.,  F(0)  is  changed  by 
Powell’s  method  so  that  conditions  (4.7)  are  satisfied  to  a  required  accuracy.  It  was  found  convenient 
in  this  example  to  represent  F (0)  by  a  cosine  Fourier  series  E™0Fj  cos  i0,  say,  where  m  is  1  or  2  for 
the  circular  cone  at  small  angles  of  incidence;  by  this  representation  the  number  of  unknowns  is 
reduced  and  so  the  work  involved  in  finding  the  partial  derivatives  in  (2.11)  is  reduced  and  also  Powell’s 
method  is  more  efficient. 

It  is  important  in  this  example  to  have  a  good  estimate  of  the  shock  shape  F(0).  To  be 
always  sure  of  a  good  estimate,  a  situation  is  first  considered  for  that  of  a  circular  cone  r  =  G(0)  =  const 
which  is  at  incidence  a  =  0  deg.  For  this  case  the  flow  is  axisymmetric  and  the  problem  is  easily  solved. 
A  situation  is  next  considered  which  has  a  small  perturbation  from  the  circular  cone  at  zero  incidence 
(either  a  small  change  to  the  body  shape  G(0)  or  a  small  change  in  incidence  may  be  considered] .  In 
this  case  the  estimate  for  F(0)  is  taken  to  be  that  obtained  for  the  first  case  of  the  circular  cone  at  zero 
incidence.  The  solution  for  this  small  perturbation  is  then  found  by  the  method  of  lines.  Next  a  larger 
perturbation  of  body  shape  or  incidence,  which  is  proportional  to  the  first  perturbation,  is  considered 
and  F(0)  is  estimated  by  extrapolation  from  the  two  previous  results.  And  so  the  technique  can  be 
continued  for  larger  proportional  perturbations  and  always  a  good  estimate  for  F(0)  is  available  by 
extrapolation  from  previous  results  at  the  smaller  perturbations.  For  example,  to  find  solutions  at 
incidence  for  the  circular  cone  whose  semi-apex  angle  is  0C,  a  solution  is  first  found  for  a  =  0  deg.,  then 
successively  for  a/0c  =  0.01,  0.1,  0.2,  0.3,  .... 

By  the  method  of  lines  it  was  possible  to  generate  solutions  for  the  circular  cone  for  relative 
incidence  a/0c  as  high  as  1.4  in  some  cases,  which  was  higher  than  relative  incidences  at  which  any 
other  theoretical  solutions  were  available.  The  only  other  methods  available  which  gave  solutions  at 
relative  incidences  greater  than  unity  (up  to  about  1.2)  are  methods  which  solved  the  full  hyperbolic 
Equations  (4.5)  [Babenko  1966,  Gonidou  1968]  and  these  methods  are  50-100  times  less  efficient 
than  the  method  of  lines. 

The  quality  of  the  results  for  a  circular  cone  can  be  seen  in  Table  VIII  which  compares 
surface  values  and  shock  shape  F(0)  on  a  circular  cone  obtained  by  MOL  with  60  =  22.5°  and  6£  =  0.1 
and  by  the  method  of  Babenko  et  al.  (1966)  with  60  =  11.25  and  5£  =  0.05,  which  solves  the  full 
hyperbolic  Equations  (4.5). 

Results  for  a  non  circular  cone  are  shown  in  Figure  4.5a  and  are  compared  with  experimen¬ 
tal  results.  The  cross  section  of  the  body  in  this  case  is  shown  in  Figure  4.5b. 

Finally  a  surface  pressure  result  of  Klunker  et  al.  (1971)  for  the  delta  wing  problem 
(compression  side)  is  given  in  Figure  4.6.  In  this  problem,  the  cross  section  of  the  wing  is  flat,  and  the 
shock  wave  is  attached  not  only  at  the  wing  apex,  but  also  along  the  leading  edges  which  are  swept 
back  50°.  The  total  velocity  is  everywhere  supersonic,  but  in  this  problem  the  conical  cross  flow  is 
also  supersonic  in  a  region  adjacent  to  the  leading  edge  (x  =  1  on  Fig.  4.6).  It  is  interesting  that  MOL 
can  be  used  without  difficulty  in  this  case,  even  though  the  conical  equations  are  of  mixed  type;  MOL 
gives  an  excellent  prediction  of  the  constant  pressure  which  occurs  in  the  hyperbolic  region  (for  the 
flat  wing)  as  well  as  in  the  elliptic  region  near  the  wing  center  line  (x  =  0  on  Fig.  4 .6),  where  the  cross  flow  is 
conically  subsonic.  In  Figure  4.6  the  MOL  results  are  compared  with  those  of  Voskresenskii  ( 1968),  who 
used  the  three-dimensional,  fully  hyperbolic,  finite-difference  method.  It  can  be  seen  that  the  two  methods 
agree  well,  and  that  remarkable  accuracy  is  obtained  by  MOL  with  only  one  intermediate  line  between  the 
wing  center  line  and  leading  edge  (i.e.,  N  =  2).  Further  MOL  results  for  delta  wing  problems  are  given 
in  Klunker  et  al.  (1971),  including  comparisons  with  experiment  and  other  calculative  methods. 
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It  may  be  noted  in  Figure  4.6  that  Klunker  et  al.  in  applying  MOL  to  the  delta  wing 
problem,  did  not  use  constant  strip  widths  but  took  more  lines  in  the  region  where  there  is  more 
variation  in  quantities,  i.e.,  near  x  =  0  in  Figure  4.6.  This  is  possible  in  their  case  since  they  approxi¬ 
mate  derivatives  by  fitting  a  fourth-order  polynomial  to  five  adjacent  points  near  to  the  point  at  which 
the  derivative  is  required. 
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FIG.  4.18  STRUCTURES  PROBLEM,  SECTION  4.2.  CONVERGENCE  AS  5 y  - 

x  =  0.5,  Sx  *  1/20,  Q  *  32. 


FIG.  4.2b  COMPARISON  OF  RESULTS  USING  N  =  3  AND  N  =  4  WITH  FINITE  ELEMENTS. 

x  =  0.5,  5x  =  1/20,  Q  =  80. 
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FIG.  4.4a  CROSS  SECTION  z  =  CONSTANT  OF  THE  FLOW  FIELD;  CONICAL  FLOW  PROBLEM 
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FIG.  4.4b 
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FIG.  4.5a  COMPARISONS  OF  EXPERIMENTAL  AND  THEORETICAL  PRESSURE 
DISTRIBUTIONS  FOR  BODY  GIVEN  BY  FOURTH  ORDER  EVEN 
COSINE  FOURIER  SERIES 


SYMMETRY 


FIG.  4.5b  BODY  SHAPE  USED  IN  CONICAL  FLOW  PROBLEM 


0  0.5  1.0 

x 


FIG.  4.6  SURFACE  PRESSURE  ON  FLAT  DELTA  WING,  M„  =  4.0,  A  =  50°, 
MOL  AND  FINITE-DIFFERENCE  HYPERBOLIC  SOLUTIONS 
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TABLE  VIII 

COMPARISONS  OF  SURFACE  VALUES  AND  SHOCK  SHAPE  BETWEEN  PRESENT 
THEORY  AND  THE  THEORY  OF  BABENKO  ET  AL.  (1966)  FOR  CIRCULAR  CONE 

=  5,  6C  =  25,  a  *  20 


0 

0 

22.5 

45 

67.5 

90 

112.5 

135 

157.5 

180 

Uj 

1.3028 

1.3167 

1.3571 

1.4205 

1.5026 

1.5956 

1.6906 

1.7741 

1.8123 

Ub 

1.3026 

1.3165 

1.3572 

1.4217 

1.5048 

1.5989 

1.6936 

1.7711 

- 

1.8129 

Vj 

0.6075 

0.6140 

0.6328 

0.6624 

0.7007 

0.7440 

0.7883 

0.8273 

0.8451 

VB 

0.6074 

0.6139 

0.6329 

0.6630 

0.7017 

0.7456 

0.7897 

0.8259 

0.8454 

W] 

0 

0.1785 

0.3464 

0.4862 

0.5890 

0.6358 

0.6205 

0.4319 

0 

WB 

0 

0.1792 

0.3446 

0.4832 

0.5818 

0.6285 

0.6075 

0.4463 

0 

Pi 

2.6838 

2.5058 

2.0418 

1.4619 

0.9277 

0.5447 

0.3172 

0.2426 

0.2522 

Pb 

2.6842 

2.5062 

2.0423 

1.4596 

0.9282 

0.5434 

0.3182 

0.2436 

0.2508 

Pi 

4.7758 

4.5474 

3.9284 

3.0946 

2.2362 

1.5288 

1.0391 

0.8580 

0.8819 

Pb 

4.7759 

4.5475 

3.9289 

3.0909 

2.2368 

1.5260 

1.0411 

0.8603 

0.8785 

Fj 

_ 

0.5919 

0.5943 

0.6029 

0.6168 

0.6395 

0.6657 

0.6960 

0.7075 

0.6920 

B 

0.5920 

0.5947 

0.6028 

0.6173 

0.6388 

0.6665 

0.6949 

0.7068 

0.6917 

Subscripts:  J  Values  obtained  by  present  method 


B  Values  obtained  by  Babenko  et  al.  (1966) 


CHAPTER  5.0  ELLIPTIC  EIGENVALUE  PROBLEMS 

5.1  Introduction 

Although  the  main  advantage  of  the  method  of  lines  has  been  in  solving  boundary  value 
problems,  particularly  the  rather  complicated  nonlinear  types  mentioned  in  Chapter  4,  we  feel  that  for 
completeness  some  coverage  should  be  given  to  eigenvalue  problems.  Thus  in  this  chapter  we  consider 
the  method  of  lines  solution  to  Helmholtz  equation  in  a  rectangle  and  compare  MOL  results  to  the 
exact  solution. 

The  instability  problems  mentioned  earlier  (Section  3.3)  are  again  encountered  in  the 
eigenvalue  case.  However  it  is  illustrated  that  accurate  results  can  be  obtained  for  the  smaller  eigen¬ 
values  using  the  usual  five  line  difference  scheme  (2.17).  We  first  consider  the  three  line  scheme  as  this 
gives  better  insight  into  the  stability  and  accuracy  expected. 

The  problem  considered  here  has  a  linear  operator  but  nonlinear  operators  would  be  handled 
in  a  similar  manner. 

5.2  The  Method  of  Lines  Applied  to  Eigenvalue  Problems 

Consider  Helmholtz’  equation  in  a  rectangle  o<x<a,  o<y<b.  Lines  are  drawn  parallel 

to  the  x  axis  and  numbered  0,1, . N.  Because  of  symmetry  we  consider  only  the  region 

o  <  x  <  a/2,  o  <  y  <  b/2.  Lines  o,  N  correspond  to  y  =  0  and  y  =  b/2  respectively.  Writing  h  =  b/2N 
and  replacing  32t///3y2  on  the  nth  line  by  the  three  line  finite  difference  scheme  (2.5),  for  example,  we 
have 


dpn  « Pn+1  -  24>n  +  <//„_, 

—  +  -  +  k2i//n  =  0  (5.1a) 

dx  h2 


where 


d'/'n 

-—  =  pn  n  =  1,2  ...  N  (5.1b) 

dx 


We  apply  symmetry  at  y  =  b/2  by  making  <An+i  =  i  •  The  initial  conditions  at  x  =  0  are 


\pn  =  0  n  =  1,2  .  . .  N 


We  can  integrate  (5.1)  from  x  =  0  to  x  =  a/2  provided  we  have  values  for  k  and  for  the 
gradients  pn  for  n  =  2,3  .  .  .  N.  The  gradient  P]  is  kept  fixed  throughout  since  we  want  to  avoid  the 
trivial  solution  =  0.  Note  that  this  fixing  of  a  gradient  +  0  may  prevent  finding  the  eigenvalue 
corresponding  to  an  eigenvector  which  has  zero  gradient  at  that  point.  For  example  if  we  fixed 
Pn_|  f  0  when  N  =  3  then  we  could  not  calculate  the  eigenvalue  corresponding  to 


mTrx  37ry 

\p  =  A  sin - sin - 

a  b 
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since 


di^„(x  =  0) 
dx 


37ry„ 

sin - 

b 


Anur  3jt(N  -  l)h 

- sin - 

a  Nh 


=  0 


when  N  =  3.  By  increasing  N  to  4,5  .....  though,  we  could  calculate  the  eigenvalue.  Choosing  P]  f  0 
for  this  particular  problem  does  not  lead  to  any  eigenvalue  loss. 


Since  the  problem  is  linear  for  a  fixed  k  we  proceed  as  follows.  Fix  k  in  an  outer  loop.  In 

N  [”  N 

the  inner  loop  calculate  p2  .  . .  .  pN  so  that  the  sum  of  squares  of  residuals  Ir?  =  Ip- 

i  =  l  i=l 

is  minimized.  The  generalized  least  squares  minimization  leads  to 


at  x  =  a/2 


N 

2 

J=2 


3rk  3rk| 


6F: 


N 


aiv 


i  =  2 . N 


(5.2) 


where  Fj  =  pj  (x  =  0)  and  5Fj  is  the  change  to  be  made  in  Fj  so  that  Sr?  is  minimized.  Since  rk  is  linear 
with  respect  to  (F2  .  .  .  .  FN ),  (5.2)  is  exact  and  the  partial  derivatives  can  be  found  exactly  (within 
round  off  errors)  by  differences 


3*k  _  rk(F2  •  .  •  ,  ^  +  AFj,  •  .  .  FN)  -  rk(F2  .  .  .  ,  Fj, ,  ,  .  FN) 
__  _  —  __  __ 


(5.3) 


where  AFj  is  theoretically  any  value  but  the  authors  usually  select  lCT6Fj  if  Fj  f  0  since  a  large  dis¬ 
turbance  to  F;  may  prevent  integration  to  x  =  a/2  to  find  rk.  Also  if  the  operator  is  non  linear  then 
(5.3)  is  a  good  approximation  only  if  AFj  is  small.  The  first  estimate  of  (Fj)  is  not  too  critical  since 
the  above  method  theoretically  yields  the  exact  minimum  in  the  inner  loop.  However  we  have  to 
select  values  which  enable  us  to  integrate  to  x  =  a/2  without  the  solution  blowing  up.  In  the  example 
of  this  chapter  we  consistently  set  F,  =  F2  =  . .  . .  FN  =  1  at  the  start  of  the  inner  loop,  used 
AFj  =  10“®Fj,  and  did  not  encounter  any  difficulties.  Fj  was  then  kept  fixed  on  1  and  integrations 
made  using  the  small  perturbations  to  F2  .  .  .  .  FN  in  turn.  Having  found  the  partial  derivatives  accord¬ 
ing  to  (5.3)  we  then  substitute  into  (5.2)  and  by  Gaussian  elimination  find  8Fj. 

In  the  outer,  k,  loop  we  can  either  select  values  for  k2  in  some  consistent  manner. 


k2  =  c(d)e  (5.4) 


say,  or  we  can  carry  out  a  one  dimensional  minimization.  The  method  used  by  the  authors  is  a  com¬ 
bination  of  both  these  possibilities.  First  c  and  e  are  selected  (for  example  we  chose  c  =  0,  e  =  65  in 
the  square  membrane  problem  to  follow)  and  the  step  d  is  selected  in  a  manner  depending  on  the 
problem  at  hand  (for  the  square  membrane  we  chose  d  =  1  since  we  knew  there  were  many  roots 
expected  in  the  range  o  <  k2  <  65).  For  a  particular  N  we  then  plot  R  (=  (2r?)Vi)  against  k2.  The 
function  R  will  have  a  cusp  at  each  root  (Fig.  5.1  for  example)  thus  giving  approximations  kf ,  k2 
to  the  eigenvalues.  We  may  then,  if  required,  calculate  the  roots  more  accurately  by  minimizing  Sr2 
with  respect  to  X  where  X  is  given  by 


■  ■■  '  ?' 


'V, 
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k2  =  '/4k,2  +  Kj2)  +  'h(k*  -  kj2)  sin  X 


(5.5) 


and  k2  is  the  ith  eigenvalue.  k(2  and  k(2  are  lower  and  upper  limits  on  k2  chosen  from  the  graphical 
inspection  above.  The  authors  preferred  this  semi-automatic  method  but  a  completely  automated 
method  could  be  devised  in  which  the  computer  would  locate  each  root  approximately  and  then  refine 
the  approximation. 

If  we  suspect  the  presence  of  multiple  roots  which  numerically  are  roots  which  may  be  very 
close  but  not  coincident  we  can  then  minimize 


Sr2/(k2  -  kf2)2  (5.6) 


in  the  ranges  kr  <  k2  <  k^  and  k^<  k2  <  k;2  where  k(2  is  the  eigenvalue  already  found.  Triple  roots 
can  be  found  in  a  similar  manner. 

The  method  has  been  described  above  using  a  three  point  approximation  to  the  derivative 
i//yy.  But  in  order  to  use  as  few  lines  as  possible  it  is  better  to  use  finite  difference  approximations 
having  greater  accuracy.  Application  in  this  case  is  exactly  the  same  as  outlined  above. 

We  next  apply  the  method  to  a  square  membrane  using  the  three  point  scheme  (for  which  an 
analytical  solution  is  available)  and  then  use  a  five  point  scheme  for  which  a  semi-analytical  solution  is 
available. 


5.3  Helmholtz’  Equation  in  a  Square 

To  investigate  the  accuracy,  stability  and  convergence  of  MOL  we  consider  a  simple  example 
with  a  known  solution  and  solve  the  MOL  ordinary  differential  equations  first  analytically  and  later 
numerically. 


5.3.1  The  Three  Line  Scheme 

We  consider  a  square  of  side  n  in  which  we  want  to  solve  (v2  +  k2)\p  =  0  subject  to  \p  =  0 
on  the  boundary.  The  MOL  representation  of  Helmholtz’  equation  using  the  three  line  scheme  can  be 
written 


*n+l  "  2tn  +  '/'n-l 

+  -  +  k2*n  =  0  (n  =  1,2  ...  N)  (5.7) 

h2 


where  h  =  jr/2N  and  lines  are  considered  parallel  to  the  x  axis.  Lines  o,  N  are  equivalent  to  y  =  0  and 
y  =  ir/2  respectively.  The  primes  denote  differentiation  with  respect  to  x.  We  consider  only  the  region 
o  <  x  <  7t/2,  o  <  y  <  tt/2  because  of  symmetry  (or  antisymmetry). 


The  boundary  conditions  are 


x  -  0:  \Pn  =  0  (n  =  1,2  .  . .  N) 


(5.8a) 
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x 


jr 

2' 


4/'n  =  0  for  symmetry 


(n  =  1,2  ...  N) 


[or  »,f/n  =  0  for  antisymmetry] 


(5.8b) 


i/'n  +  i  =  tf'N-i  for  symmetry 


[or  +  \  for  antisymmetry] 


(5.8e) 


^o=0 


(5.8d) 


It  can  be  shown  that  the  general  solution  to  (5.7)  subject  to  the  boundary  conditions  (5.8c) 
and  (5.8d)  is 


4  mmr 

K  =  2  sin  — —  [Am  sin  ^mx  +  Bm  cos /imx]  (5.9) 

m=l  ,3,5.  .  .  2N 


where  =  k2  -  (4/h2)  sin2  m7r/4N.  Since  the  solution  is  real  it  follows  that  Am  is  imaginary  if 
k  <  (2/h)  sin  rmr/4N.  We  now  apply  the  remaining  boundary  conditions  (5.8a)  and  (5.8b)  giving 


nm7r 

2  sin  — —  Bm  =  0  (n  =  1,2  . .  .  N)  (5.10a) 

m  2N 


and 


2 

m 


sin 


nm7T 

2N~ 


Am  cos  — 


®m  sin  “ 


0  (n  =  1,2...  N) 


(5.10b) 


The  system  (5.10a),  assuming  that  det  (anC)  f  0  where  an6  =  sin  n(28-l)7r/2N,  shows  that  B]  =  B3 
=  .  .  .  =  B2N_]  =  0  and  the  system  (5.10b)  becomes 


2  sin 
m 


nmir 

2N 


Mm* 

cos  —  nmAm 


=  0 


(5.10c) 


For  a  non  trivial  solution  we  require  HmAm  f  0  for  some  m  (see  Equation  (5.9))  and  hence  require  the 
determinant  of  system  (5.10c)  to  be  zero  i.e. 
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2N-1 

n  cos  nm 

m=l ,3,5  .  .  . 


7T 

2 


(det  anli )  =  0 


or  since  we  assume  de^a,^)  f  0 


and  so 


cos  Mm  -  =  0  for  m  = 

=  p  where  p  is  one  of  the  integers  1,3,5.... 


1,3 _ or  2N-1 


Thus 


p  -  1,3,5  .  .  . 
m  =  1,3, .  . .  2N-1 


(5.11) 


(5.12) 


giving  the  MOL  analytical  approximation  to  the  exact  eigenvalues  which  in  the  symmetrical  case  are 


k2  =  p2  +  m2 


P 

m 


1,3,5  .  . . 
1,3,5  .  .  . 


Since  h  =  tt/2N  the  difference  is 


16N2 


ITT 


mV 


16N2 


__4  4 
m  ir 

768N4 


_  mV 

48N2 


+  0(N-4) 


(5.13) 


(5.14) 


indicating  accuracy  of  0(h2)  as  expected.  It  can  be  seen  from  (5.14)  that  for  the  smallest  eigenvalue 
(m  =  1)  the  error  will  be  <  10~3  if  N  >  15.  For  m  =  3  using  the  same  number  of  lines  the  error  would 
increase  by  a  factor  81  and  higher  eigenvalues  would  rapidly  increase  in  error.  Values  of  k2-p2  found 
from  the  three  point  scheme,  formula  (5.12),  are  listed  in  Table  IXa  for  N  =  3,4, 5, 6, 7  and  8.  This 
table  illustrates  that  accuracy  for  m  =  1  is  fairly  good  when  N  =  8  but  that  higher  eigenvalues  are 
poorly  approximated. 

The  conclusion  of  this  subsection  is  that  many  lines  have  to  be  used  with  the  three  point 
scheme  in  order  to  obtain  a  reasonable  accuracy  even  for  the  smallest  eigenvalue.  Using  such  a  large 
number  of  lines  would  undoubtedly  lead  to  a  large  instability  when  applying  MOL  numerically  and  so 
we  seek  a  difference  scheme  which  is  more  accurate  than  the  three  line  scheme.  This  will  enable  us  to 
use  fewer  lines  and  so  help  to  minimize  the  effect  of  the  instability.  In  the  next  section  we  next 
consider  such  a  scheme. 

Note  that  the  eigenvector  for  the  three  point  scheme  is 


*n 


,  nmir 

A  sin  px  sin - 

2N 


A  sin  px  sin  myn 


which  is  exact  in  this  case. 
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5.3.2  The  Five  Line  Scheme 

In  this  case  i//yy  in  Helmholtz’  equation  is  replaced  by  formula  (2.17b)  on  interior  lines  and 
by  formula  (2.20)  on  the  line  adjacent  to  y  =  0.  As  with  the  three  line  scheme  the  resulting  problem  of 
finding  the  eigenvalues  can  be  reduced  to  finding  the  eigenvalues  of  a  matrix.  However  these  cannot  be 
found  in  closed  form  as  before  but  are  found  numerically. 

Applying  the  five  line  scheme  to  Helmholtz’  equation  results  in  a  set  of  ordinary  differential 
equations  of  second  order  in  x  written  at  each  line  1,2  ....  N.  On  seeking  a  solution  of  the  type 
t^n  -  an  sin  px  we  have 


a2  _  ai  i  75 

(k2  -  p2)ai - —  +  "j-(a3  -  a,)  -  !4(a4  -  a,) 


0.125  1.25 

+  Z  <a5  -  ai )  +  ~T~  («o  -  ai )  =  0 


h2(k2  -  p2)a„  +  ^(an+1  +an_,  -  2an)  -  ^a„+2  +  a„_2  -  2a„)  =  0 


(5.15a) 


n  =  2,3  ...  N 


(5.15b) 


with  boundary  values  aQ  =  0,  aN+1  =  aN  _  1  and  aN+2  =  aN  _2  for  symmetry.  To  give  a  non  trivial 
solution  the  determinant  of  the  matrix  of  system  (5.15)  should  be  zero.  For  example  using  N  =  3  we 
require 


7 

T-  - 
6 


31 

T - 

12 


5 


where  T  =h2(k2  -  p2)  i.e.  we  require  the  eigenvalues  of  the  above  matrix  written  without  T  on  the 
main  diagonal.  If  X( ,  X2  . . . .  XN  are  the  resulting  eigenvalues  then 


k2  -  p2  =  -  —  i  =  1 ,2  .  .  .  N 
h2 


(5.16) 


Now  p  =  1,3,5  ....  to  satisfy  the  boundary  conditions  i^n  =  0  at  x  =  0  and  \pn  symmetrical  about 
x  =  n/2.  Thus  the  eigenvalues  k2  are  determined  from  (5.16).  Values  of  k2-p2  from  Equation  (5.16) 
are  listed  in  Table  IXb  for  N  =  3, 4,5, 6, 7  and  8.  This  table  illustrates  the  accuracy  to  be  expected  from 
the  five  line  scheme  since  the  exact  k2  -p2  should  be  equal  to  m2  for  m  =  1,3, ....  2N-1.  The  accuracy 
of  the  eigenvalues  corresponding  to  m  =  1  and  3  is  very  good  while  the  m  =  5  and  7  eigenvalues  have 
about  1  %  accuracy  at  N  =  8.  Higher  eigenvalues  are  poorly  predicted.  Note  that  convergence  is  rapid, 
for  example  the  errors  for  the  m  =  5  eigenvalue  are  about  30%,  5%,  3%,  1.5%  and  0.8%  with  N  =  4, 5, 6,7 
and  8  respectively.  This  feature  is  to  be  expected  since  the  error  is  0(N“4). 

To  compare  the  above  MOL  results  with  those  obtained  by  finite  differences  we  can  use  the 
formula  given  in  Isaacson  and  Keller  (1966)  for  the  finite  difference  solution  of  accuracy  0(h2)  using 
an  N  X  N  mesh  on  a  quarter  of  the  square  of  side  n  i.e. 


k2 


P»  .  2 

—  +  sin^ 
4N 


c.f.  Equation  (5.12).  For  example,  to  achieve  1%  accuracy  for  the  eigenvalue  corresponding  to  m  =  5 
and  p  =  1  would  require  a  22  X  22  mesh. 


5.3.3  Numerical  Solution  by  the  Five  Line  Scheme 

The  three  line  and  five  line  schemes,  solved  analytically  above,  illustrate  the  accuracy  one 
might  expect  to  obtain  using  MOL.  We  now  proceed  to  the  numerical  procedure  for  solving  the 
equations  and  use  the  five  line  scheme.  But  first  we  inspect  the  numerical  instability  which  is  similar 
to  the  instability  encountered  in  boundary  value  problems  (see  Section  3.3). 

As  can  be  seen  from  the  three  line  analytic  solution  (5.9,  5.10)  we  should  have  Bin  =  0  and 
Am  =  0  except  for  the  one  A,  say  a,  corresponding  to  nm  =  p  i.e. 


,  nm7r 

'Pn  =  a  sin  sin  px 


is  the  exact  solution  for  three  point  difference  scheme.  However  numerically,  since  we  treat  the 
problem  as  an  initial  value  problem,  we  will  not  have  exactly  Bm  =  0  and  Am  =  0.  Instead  these  will 
take  on  small  values  perhaps  of  order  10"10.  And  so,  unwanted  terms  in  expression  (5.9)  for  \pn  will 
be  present  and  cause  the  inherent  instability.  This  is  likely  to  be  significant  if  k  <  (4N/jt)  sin  m7r/4N 
since  then  there  is  a  term  present  in  (5.9)  containing  the  factor 


which  for  m  =  2N-1  is 
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and  when  seeking  the  eigenvalue  corresponding  to  m  =  1  for  example,  with  N  =  8  and  p  =  1  say,  a 
factor  can  be  produced  of  size  sinh  lOx.  At  x  =  tt/2  therefore,  even  though  A2N_ ,  may  be  small,  the 
solution  could  be  affected  by  this  factor.  Clearly  then  the  number  of  lines  we  use  cannot  be  too  large 
and  it  will  be  seen  later  that  the  N  =  8  results  are  slightly  affected  by  the  instability  while  N  =  7  results 
appear  unaffected. 

In  the  previous  analytical  solution  of  the  MOL  ordinary  differential  equations  we  obtained  a 

solution  explicitly  involving  x  in  the  form  sin  px  (p  =  1,3 . ).  Our  choice  of  integration  step  length 

5x  will  now  determine  how  accurately  this  factor  is  determined.  Suppose  we  require  the  error  in  sin  px 
to  be  less  than  10  3  at  x  =  n/2.  The  fourth  order  Runge  Kutta  and  Hamming’s  predictor  modifier 
corrector  (Hamming,  1959)  methods  used  in  the  computations  yield  errors 


E  < 


5  xs 
120 


fv(£) 


on  each  step  where  x  <  £  <  x  +  6x  and  f  is  the  function  of  x  being  approximated.  Integration  from 
0  to  7r/2  in  n  steps  therefore  gives  a  total  error 


from  which  follows  n  >  3p5/4  or,  for  p  =  5,  n  >  24.  Since  this  is  an  upper  bound  on  the  error  it  was 
decided  to  use  n  =  20  steps  from  0  to  tt/2  knowing  then  that  this  would  produce  sufficiently  accurate 
results  for  p  =  1,3  or  5. 

It  was  decided  to  seek  eigenvalues  in  the  range  0  <  k2  <  65  using  N  =  3,4,5, ...  as  high  as 
possible  before  inherent  instability  swamped  the  solution;  this  turned  out  to  be  N  =  9.  Since  we 
expected  a  large  number  of  eigenvalues  in  the  above  range  and  also  some  double  and  triple  roots  we 
selected  an  increment  6k2  =  1  to  locate  roughly  the  eigenvalues  by  the  process  mentioned  earlier,  see 
Equation  (5.4).  The  resulting  plots  of  (Sr2)w  against  k2  are  shown  in  Figure  5.1  where  r(  is  the  residual 
in  \p\  at  x  =  n/2  after  carrying  out  the  minimization  (5.2).  The  computation  time  for  N  =  3,4  ....  8 
inclusive  was  2  mins  on  an  IBM  3032  computer. 

The  figures  show  that  some  roots  are  repeating  for  all  values  of  N  and  therefore  are  reason¬ 
ably  accurate.  Other  roots  however  are  not  repeating  and  these  are  poor  approximations  to  eigenvalues. 
For  example,  with  N  =  3,  m  =  5  and  p  =  1  we  should  have  k2  =  p2  +  m2  =  26  but  obtain  a  root  k2  — 14.3 
which  is  consistent  with  values  obtained  in  Table  IXb.  Having  obtained  these  approximations  we  then 
locate  the  roots  more  accurately  by  minimization  using  Equation  (5.5). 

If  we  suspect  double  or  triple  roots,  for  example  near  k2  =  50,  we  find  one  root  in  the  usual 
manner  and  then  minimize  the  function  shown  in  (5.6). 

Table  IXc  lists  the  values  of  km2  corresponding  to  (m,p)  =  (1,1)  (3,3)  (1,7)  (7,1)  and  (5,5) 
obtained  by  this  technique  for  N  =  4,6,7  and  8.  It  can  be  seen  that  the  values  are  consistent  with  those 
of  Table  IXb  (add  on  p2  to  those  values  listed  where  p  =  1,3,5, 7). 

The  N  =  8  solutions  are  affected  by  the  inherent  instability.  This  is  demonstrated  by  Sr2 
being  as  large  as  0.0015  when  k2  =  2.0032.  The  next  eigenvalue  k32  is  better  than  that  predicted  in 
Table  IXb  and  therefore  is  likely  to  be  a  fortunate  effect  of  the  instability.  The  residuals  in  the  (3,3) 
mode  were  much  smaller  than  those  of  the  (1,1)  mode;  this  is  caused  by  greater  instability  in  the  latter 
case. 
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The  accuracy  of  the  results  can  be  seen  to  be  very  good  for  the  (1,1)  mode  while  for  the 
higher  modes  the  maximuni  error,  excluding  m  =  7,  is  about  0.8%.  The  modes  corresponding  to  m  =  7 
are  not  listed  for  N  <  8  since  they  clearly  are  not  very  accurate  even  when  N  =  7  (Table  IXb).  The 
(1,7)  mode  is  not  as  accurate  as  the  (1,1)  mode  —  greater  accuracy  for  the  (1,7)  mode  would  require  a 
smaller  fix. 

The  local  minimizations  used  to  produce  Table  IXc  took  about  one  min.  on  an  IBM  3032. 
Powell’s  (1964)  one  dimensional  minimization  was  used. 
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TABLE  IXa 

EIGENVALUES  k2-p2  FOR  THREE  LINE  SCHEME  GIVEN  BY  FORMULA  (5.12) 


m  =  1 

3 

5 

7 

N  =  3 

0.9774 

7.30 

13.61 

4 

0.9872 

8.01 

17.93 

24.95 

5 

0.9918 

8.35 

20.26 

32.18 

6 

0.9943 

8.55 

21.63 

36.73 

7 

0.9958 

8.67 

22.48 

39.72 

8 

0.9968 

8.74 

23.06 

41.76 

EXACT 

1 

9 

25 

49 

TABLE IXb 


EIGENVALUES  k2  p2  FOR  FIVE  LINE  SCHEME  GIVEN  IN  FORMULA  (5.16) 


m  -  1 

3 

5 

7 

9 

11 

N  =  3 

BIBf 

8.511 

13.28 

4 

9.461 

16.74 

30.1 

5 

1.0000 

9.150 

23.64 

30.6 

50.4 

6 

1.0000 

9.046 

25.66 

37.3 

54.1 

75 

7 

1.0000 

9.014 

25.42 

46.3 

56.5 

83 

8 

1.0000 

9.003 

25.21 

49.34 

65.9 

86 

EXACT 

1 

9 

25 

49 

81 

121 

TABLE  IXc 

NUMERICAL  EVALUATION  OF  THE  EIGENVALUES  k2  p  OF  HELMHOLTZ’ 
EQUATION  IN  A  SQUARE  OF  SIDE  i r 


k2 

Kll 

k2 

k33 

k2 

KI7 

k2 

*71 

k2 

*55 

N  =  4 

2.0001 

18.461 

49.934 

6 

2.0000 

18.046 

49.934 

7 

2.0000 

. 

49.933 

50.41 

8 

2.0032 

49.940 

50.34 

50.20 

EXACT 

2 

18 

50 

50 

50 
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(2r 


FIG.  5.1  RESIDUALS  (Erf  I*  v«.  k2  FOR  SQUARE  OF  SIOE  *. 
(m,p)  ON  k2  AXIS  REFERS  TO  EIGENVALUE  k2,  FROM  TABLE  IXb 
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APPENDIX  A 


The  Method  of  Integral  Relations 


As  mentioned  in  the  Introduction,  there  is  another  semidiscrete  method  which  has  been 
widely  used  in  aerodynamic  problems,  called  the  method  of  integral  relations  (MIR).  Holt  (1977) 
reviews  the  aerodynamic  applications  of  the  method,  but  theoretical  treatments  or  studies  of  asymptotic 
convergence  are  rare. 

The  method  is  usually  applied  to  systems  of  first-order  partial  differential  equations.  As  in 
MOL,  the  region  is  considered  to  be  divided  into  strips  which  are  parallel  to  one  co-ordinate,  x  say.  The 
equations  are  partially  integrated  with  respect  to  the  other  co-ordinate,  y,  to  obtain  an  approximate 
system  of  ordinary  differential  equations.  The  partial  integration  is  performed  explicitly  by  assuming 
an  appropriate  y  dependence  of  the  integrands.  Most  applications  have  used  for  this  purpose  a  polyno¬ 
mial  whose  degree  increases  proportionally  to  the  number  of  strips.  The  algebraic  development 
required  for  this  procedure  becomes  very  cumbersome  for  N  >  2,  so  several  investigators  have  used  a 
linear  y  dependence  from  one  line  to  the  next  in  order  to  obtain  a  simple  recursive  form  for  a  system  of 
ordinary  differential  equations.  We  will  consider  the  application  of  the  latter  procedure  to  the  example 
of  Section  3.3,  since  an  explicit  solution  can  then  be  found. 

First  the  substitution  ipx  =  u,  i£y  =  -  v  is  made  in  Equation  (3.4)  to  obtain  the  Cauchy- 
Riemann  equations 


ux  -  vy  =  0;  vx  +  uy  =  0. 


With  the  polygon  approximation  for  the  y  dependence,  the  partial  integration  with  respect  to  y  yields 


vn+1  -  v„ 

l/2(un+1  +  un) - - -  *  0, 


un+i  "  un 

^(v„  +  vn+1)  +  - - -  =  0, 


where  the  notation  is  similar  to  that  introduced  in  Equations  (3.8)  and  (3.11).  Differentiation  of  (A2) 
with  respect  to  x  and  manipulation  yields: 


‘/4(u”+l  +  2u'n  +  u;_t)  +  —  (u„+l  -  2un  +  u,,.,)  =  0, 

h2 


producing  a  tridiagonal  system  for  the  x  derivatives. 

The  solution  of  the  above  system  is  in  the  identical  form  of  Equations  (3.20)  and  (3.22), 
except  that 


\Hm\  =  (2N/b)tan(mff/4N), 


t 
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cosh  0  - 


o  -  (i  ♦  ^/(i  - 

\  4N2//\  4N2/ 


From  Equation  (A5)  we  see  that  the  largest  eigenvalue  is 


M2N-1  38  8N2/irb, 


giving  Equation  (3.25)  for  the  instability  factor  in  =  un(‘/2).  Expansion  of  the  hyperbolic  cosine 

function  as  in  Section  3.4  shows  that  the  MIR  discretization  error  is  0(N-2);  but  MIR  is  clearly  inferior 
to  MOL  from  the  viewpoint  of  the  size  of  the  eigenvalues  pm.  That  is,  the  MOL  eigenvalues  grow 
linearly  with  N,  while  the  MIR  eigenvalues  are  quadratic  in  N.  Further,  the  extra  complication  of  the 
tridiagonal  system  for  the  x  derivatives  has  been  added  without  gaining  the  benefit  of  decreased 
y-truncation  error,  contrary  to  the  scheme  of  Appendix  B. 
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APPENDIX  B 


Tridiagonal  MOL  Systems  With  Accuracy  0(N~4) 


As  alternatives  to  the  five-point  difference  schemes  (2.17)  we  present  here  two  schemes  with 
the  same  accuracy  0(N-4)  which  involve  only  three  adjacent  lines.  The  schemes  are  not  general  but  can 
be  derived  in  certain  cases  as  follows. 


Poisson  Equation 

We  have  in  the  notation  of  this  paper  from  Taylor-series  expansions 


320 n  h4  340n 

*Pn+t  ~  20n  +  0n_,  =  h2  - —  +  — - —  +  0(h6) 

3y2  12  3y4 


and  also 


<*2K+\ 

by2 


b2\pn  a2^„_1 

2 -  +  - 

by1  9y2 


+  0(h4). 


(Bl) 


(B2) 


On  eliminating  340n/3y4  from  (Bl)  and  (B2)  and  substituting 


a2  020k 

— —  = - +  f(x,yk)  (k  =  n  -  1,  n,  n  +  1) 

3y2  3x2 


from  Poisson’s  equation,  we  obtain 


^(Cl  +  10K  +  K-l)  +  h  2(^n+l  "  20n  +  0„_ t) 


=  —  (f(x,  yn+i)  +  10f(x,yn)  +  f(x,  y„_i)), 


giving  a  tridiagonal  system  of  equations  with  accuracy  0(h4)  i.e.,  0(N  4). 


(B3) 


(B4) 


In  fact  the  system  (B4)  can  be  solved  analytically  in  certain  cases  and  in  particular  we  refer 
to  the  example  of  Section  3.3  and  obtain  the  solution.  It  can  be  shown  that  this  solution  is  identical 
to  (3.20)-(3.22)  except  that 
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2N  rmr 
Mm  =  —  sin  — 
m  b  4N 


/  1  .  2  mirV 

1 1  -  —  sm:  —  ) 
\  3  4N  / 


cosh  6 


('  •  if )/(' 


With  this  scheme  the  y-truncation  error  is  reduced  by  two  orders  of  magnitude  compared  to  (2.5) 
while  the  largest  eigenvalue  M2N-1  *s  increased  only  by  a  constant  factor  of  about ^/(l. 5).  We  may 
expect  to  obtain  results  of  accuracy  comparable  to  the  five-point  scheme  (2.17). 


First-Order  Equations 


Consider  the  first-order  equations 


3Pi  3Q; 

—  +  —  +  Ri  =  ° 

3x  3y 


i  =  1,  2, . .  .  ,  m,  where  Pj,  Qj  and  R;  are  linear  or  nonlinear  functions  of  the  independent  and  depend¬ 
ent  variables,  e.g.,  Pj  =  Pj(x,  y,  u, ,  u2, .  .  . ,  um).  For  instance  the  governing  equations  for  two- 
dimensional  flow  can  be  written  in  the  above  form. 

Dropping  the  i  subscript  and  using  Qn  to  denote  the  value  of  Q  on  the  n-th  line  for  any  one 
of  the  i  values  in  Equation  (B7)  we  have  from  Taylor-series  expansions 


Qn+l  "  Qn-l  dQn  h2 


3y  6 


+  0(h4) 


and  also 


9Qn+l  ^Qn  9Qn-l  d3Qn 

— - 2  —  +  —  =  h  — -  +  0(h4). 

3y  3y  3y  9y3 


Eliminating  33Qn/3y3  from  (B8)  and  (B9)  and  substituting 


(3Qk/3y)  =  -  Rk  -  (3Pk/3x)  (k  =  n-  l,n,  n  +  1) 


(BIO) 


J 
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i(p;+J  +  4p;  +  p;_,)  + 


Qn+l  Qn-1  1 

- - -  +  -(Rn+1  +  4Rn  +  Rn_t)  =  0,  (Bll) 

zn  d 


where  P, 
order  h 


?'v  =  dPk  /dx.  Hence  a  tridiagonal  system  for  the  x  derivatives  is  obtained  which  has  error  of 
\  i-e.,  0(N-4). 


